利用Python获取地理信息并可视化

Tech

网络中有丰富的地理信息数据资源，此处介绍我采用Python工具实现地形地貌、行政区划、城市交通路网等数据获取及可视化的部分实例，供大家参考学习。

1. 简易地形底图绘制

1.1 WMS数据库

针对较大范围的研究区域，我采用WMS数据库的地形底图，依四角坐标下载特定尺寸的图片，并展示。此处以东亚地区为例，逐步展示。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# 1. Set the left, bottom, right, top locations
llx,lly, urx,ury = [104.211,29.150,105.523,30.150]

# 2. Download the original file
target = 'http://www2.demis.nl/wms/wms.asp?Service=WMS&WMS=worldmap&Version=1.1.0&Request=GetMap&BBox=%s,%s,%s,%s'%(llx,lly, urx,ury) +\
         '&SRS=EPSG:4326&Width=800&Height=600&Layers=Bathymetry,Countries,Topography,Ocean%20features,Hillshading,Borders,Waterbodies,Coastlines&Format=image/gif'
path_name = './Terrain.gif'
urllib.urlretrieve(target, path_name) 
	
# 3. Plot the terrain map of East Asia
from scipy.misc import imread
fig = plt.figure(figsize = (8,6))
ax = plt.subplot()
m = Basemap(projection='cyl', resolution='l',llcrnrlon = llx,llcrnrlat=lly,urcrnrlon = urx,urcrnrlat=ury)
img = imread(path_name)
plt.imshow(img, zorder=1,extent=[llx,urx,lly,ury ] )
m.drawparallels(np.arange(-90., 120., 10.), labels=[1, 0, 0, 0],fontsize =18,linewidth= 0.02)
m.drawmeridians(np.arange(45,180,10), labels=[0, 0, 0, 1],fontsize = 18,linewidth= 0.02)
plt.subplots_adjust(left=0.015, bottom  = 0.05, top = 0.99,right=0.99, wspace = 0.1,hspace=None)

出图的样例，如下所示。

enter image description here

1.2 SRTM geotiff数据

对于较小的地理环境，直接采用WMS格式图片清晰度不够。此处介绍下载并绘制原始DEM地形文件的方法。
通过elevation package进行tif格式地形文件的下载

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
# 以四川某市及周边地区为例
## 1. Download the geotiff file
!pip install elevation  
!eio --product SRTM3 clip -o DEM.tif --bounds 104.211 29.150 105.523 30.150  # left bottom right top locations
	
## 2. Read the tif file as 2-d array
pathToRaster = r'./DEM.tif'
raster = gdal.Open(pathToRaster,  gdal.GA_ReadOnly)
dem = raster.GetRasterBand(1).ReadAsArray()
dem = dem[::-1]
	
	
## 3. Generate a terrain coloarmap 
import matplotlib.colors as colors
def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100):
	new_cmap = colors.LinearSegmentedColormap.from_list(
       'trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=minval, b=maxval),
		cmap(np.linspace(minval, maxval, n)))
	return new_cmap
cmap = plt.get_cmap('terrain')
terrain_cmap = truncate_colormap(cmap, 0.2, 0.7)
	
## 4. Plot it!
llx,lly, urx,ury = [104.211,29.150,105.523,30.150]
fig = plt.figure(figsize = (8,5))
ax = plt.subplot()
m = Basemap(projection='cyl', 	resolution='l',llcrnrlon = llx,llcrnrlat=lly,urcrnrlon = urx,urcrnrlat=ury)
	dem_plot = map.imshow(dem, extent=[llx,urx,lly,ury ], cmap =plt.get_cmap('terrain'),vmin = 0 ,vmax = 900)# terrain_cmap
map.drawparallels(np.arange(-90.0, 120., 0.5), labels=[1, 0, 0, 0],fontsize = 12,linewidth = 0.1,color = '#FFFFFF', zorder = 10)
map.drawmeridians(np.arange(-180.0, 180., 0.5), labels=[0, 0, 0, 1],fontsize = 12,linewidth = 0.1,color = '#FFFFFF',zorder=10)

结果如下图所示：
enter image description here

1.3 Cartopy STAMEN 图形绘制

数据源stamen是开源的地理信息服务网站，提供道路、地形等多种数据信息。在其官方网页采用手动下方式，仅能以确定中心点和缩放比例的方式进行图片下载，不能实现任意区域内图形获取。此处参照Geology and Python博客中的相关教程，以中国全域为例，展示底图获取和绘制流程。

cartopy是用以分析并可视化地理数据的Python库。
在安装时候，与系统原有的matplotlib.basemap发生冲突，原因在于二者所依赖的GEOS库不一致。在此处花费了不少时间来解决这一问题(网络上也有很多人遇到了类似的问题，如无法同时安装成功，安装后无法import以及执行cartopy有关指令，kernel会自动restart. )。

在此处，我列出了自己的安装步骤: 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
## 卸载原有的cartopy及相关库shapely，保留原有的basemap
conda uninstall cartopy    
conda uninstall shaeply
	
## 升级gdal库  
brew upgrade gdal
	
## 采用pip方式安装cartopy及shapely
pip install shapely cartopy --no-binary shapely --no-binary cartopy

## 结果仍未能成功import cartopy,原因在于shapely库未安装成功。用conda方式再安装
conda install shapely 
	
## 全部安装完毕!! d(`･∀･)b

## 开始绘图

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.feature import BORDERS
from cartopy.io.img_tiles import StamenTerrain
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

def basemap_terrain(extent):
	fig = plt.figure(figsize=(8, 6))
    ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
    ax.set_extent(extents=extent, crs=ccrs.Geodetic())
    ax.coastlines(resolution='10m')
    BORDERS.scale = '10m'
    ax.add_feature(BORDERS)
    ax.gridlines(color='.5')
    ax.set_xticks(np.arange(extent[0]+2, extent[1], 15), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(extent[2]+2, extent[3], 15), crs=ccrs.PlateCarree())
    ax.xaxis.set_major_formatter(LongitudeFormatter())
    ax.yaxis.set_major_formatter(LatitudeFormatter())
    return fig, ax

# 地图四至
x1,x2,y1,y2 = 74,134, 18,51
extent      = x1 -1, x2+1,y1-1, y2+1     
	
fig, ax = basemap_terrain(extent)
st = StamenTerrain()
ax.add_image(st, 4)

enter image description here

​
​

2. 行政区划及城市路网信息

此处以四川内江市为例，利用osmnx库，shp文件自动下载于工作文件夹内。 

1
2
3
4
5
6
7
# 行政辖区下载
import osmnx as ox
city = ox.gdf_from_place('Neijiang City,  Sichuan, China')
ox.save_gdf_shapefile(city)
city = ox.project_gdf(city)
fig, ax = ox.plot_shape(city, figsize=(3,3))

enter image description here

1
2
3
4
# 内江市东兴区公路路网下载
Road = ox.graph_from_place('Dongxing, Neijiang,  Sichuan, China', network_type='drive')
Road_projected = ox.project_graph(Road)
fig, ax = ox.plot_graph(Road_projected)

enter image description here

3. 一些有用的小指令

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
# A. 任意地点的海拔查询
!pip install geocoder
import geocoder
print geocoder.elevation ([lat,lng]).meters

# B. 两例获取某地区及周边的google地图静态图像
## B1. 卫星图
http://maps.googleapis.com/maps/api/staticmap?center=42.2341289,118.9790778&zoom=12&format=png&sensor=false&size=2400x2000&maptype=satellite&style=feature:administrative|weight:0.1|invert_lightneºss:true|visibility:off&style=feature:water|element:labels.text|visibility:off&style=feature:landscape.natural|visibility:on

## B2. 道路地图
http://maps.google.com/maps/api/staticmap?sensor=false&size=512x512&center=Brooklyn&zoom=12&style=feature:all|element:labels|visibility:off

# C. google地图下载工具
http://www.chengfolio.com/google_map_customizer#satellitemap

# D. 地形图在线下载工具
https://maps-for-free.com/
    
# E.测量两点间距离
# source: https://stackoverflow.com/questions/15736995/how-can-i-quickly-estimate-the-distance-between-two-latitude-longitude-points
from math import radians, cos, sin, asin, sqrt
def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees)
    """
    # convert decimal degrees to radians 
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a)) 
    # Radius of earth in kilometers is 6371
    km = 6371* c
    return km

#Python
中国空气质量历史数据整理(HDF5格式) MODIS卫星火点数据的简易处理与可视化分析

Kommentare

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×