前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >Python中关于底图的操作

Python中关于底图的操作

作者头像
MeteoAI
发布2019-12-25 14:35:55
3.1K2
发布2019-12-25 14:35:55
举报
文章被收录于专栏:MeteoAI

两个人,一个爱总结,一个爱技术

前言

Python基于其强大的功能越来越成为了科学利器,气象上对精细化的要求越来越高,对于底图的制作也越来越高。本人气象出身,长期用NCL画图,但是NCL对于精细化底图的支持很差(或者说因为本人不是地图学专业不明白shp文件),也不愿意学Arcgis,于是和同事小陈折腾了一系列的在Python下地图的操作。

大概有这几个部分,1)县级边界的剪切;2)添加乡镇边界;3)省市县三级边界;4)关于海洋的掩膜

先前准备

做这些之前,首先准备好CHN_adm_shp.rar文件,解压缩后有这些东西。除了知道shp文件外剩下的不懂。但是百度得知,CHN_admn0.csv类似这种的文件表示行政区的级别,比如CHN_admn1.csv到省级,CHN_admn2.csv到市级,CHN_admn3.csv到县级。在后面的代码中将利用到这些信息进行索引。具体的打开csv文件后大家结合百度,就能明白。

地址链接:

https://pan.baidu.com/s/13y493ZvJJkAQ1Omrq1dfQQ

Talk is cheap, show me your code.

代码语言:javascript
复制
from numpy import *
from netCDF4 import Dataset, date2index, num2date,MFDataset
from mpl_toolkits.basemap import Basemap,addcyclic
import os
import numpy
import numpy as np
numpy.set_printoptions( threshold=500)
import struct
from matplotlib.pyplot import *
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.colors import LinearSegmentedColormap
from matplotlib import colors
import matplotlib as mpl
import matplotlib
import matplotlib.patches as mpatches
import shapefile as shp  
def xshow(filename, nx, nz):
    f = open(filename, "rb")
    pic = np.zeros((nx, nz))
    for i in range(nx):
        for j in range(nz):
            data = f.read(4)
            elem = struct.unpack("f", data)[0]
            pic[i][j] = elem
    f.close()
    return pic
#上面的这段代码是读取dat文件用的  
fig, ax = plt.subplots(figsize=(15,15))   
m = Basemap(projection='merc',resolution='c',llcrnrlat=25.2,urcrnrlat=26.1,
                                            urcrnrlon=116.8,llcrnrlon=115.9,lon_0=120.,lat_0=90,
                                            llcrnrx=0)

landsat = xshow('20131005_ndvi_pro_jc.dat',3002,3002)
landsat1 = np.ma.masked_array(landsat, np.isnan( landsat) | (landsat<0) | (landsat>1) )
#进行县级剪裁的重要代码在下面一段,通过索引县级找出changting后,剩下的区域全部填成白色
shp_info = m.readshapefile("CHN_adm_shp\\CHN_adm3",'states',drawbounds=False,linewidth = 0.4) # 读取CHN_adm3的意义在于进行县级索引
for info, shp in zip(m.states_info, m.states):
    proid = info['NAME_3']  
    if proid == 'Changting':   
        poly = Polygon(shp,facecolor='None',edgecolor='b', lw=0.8)
        ax.add_patch(poly)
    else:
        poly = Polygon(shp,facecolor='w',edgecolor='w', lw=0.8)
        ax.add_patch(poly)
m.imshow(landsat1[::-1,:],norm=matplotlib.colors.BoundaryNorm([0,0.1,0.2,0.3 , 0.4,0.5, 0.6,0.7,  0.8, 0.9, 1.0],cm.RdYlGn.N),cmap='RdYlGn')
#这里模仿了gis的图注方式
legend_elements = [mpatches.Patch(facecolor='#a50026',label='0-0.1'),
                   mpatches.Patch(facecolor='#da362a',label='0.1-0.2'),
                   mpatches.Patch(facecolor='#f67a49',label='0.2-0.3'),
                   mpatches.Patch(facecolor='#fdbf6f',label='0.3-0.4'),
                   mpatches.Patch(facecolor='#feeda1',label='0.4-0.5'),
                   mpatches.Patch(facecolor='#ebf7a3',label='0.5-0.6'),
                   mpatches.Patch(facecolor='#b7e075',label='0.6-0.7'),
                   mpatches.Patch(facecolor='#75c465',label='0.7-0.8'),
                   mpatches.Patch(facecolor='#249d53',label='0.8-0.9'),
                   mpatches.Patch(facecolor='#006837',label='0.9-1')]
plt.legend(handles=legend_elements, loc='lower right',title='NDVI')
m.drawmapscale(116.1, 25.25, 120, 90, 30, units='km',barstyle='fancy', yoffset=0.05)
plt.show()

出图效果如下

下面我们进行第二步,在县级底图上精确到乡镇。由于上述的CHN_adm_shp.rar只能精确到县级,又下载了一份乡镇级别的shp。并在上述代码段中35行后加上这一块代码段,于是乡镇边界上去了。

代码语言:javascript
复制
#画出乡镇边界这段为这段代码
sf = shp.Reader("hetianzhen.shp",encoding="gbk")
for nn,shape in enumerate(sf.shapeRecords()):
    x = [i[0] for i in shape.shape.points[:]]
    y = [i[1] for i in shape.shape.points[:]]
x1,y1 = m(x,y)
plt.plot(x1,y1,color='black',lw=0.8)
#################################################################

出图效果如下

那么下面进行多底图加强版

代码语言:javascript
复制
from netCDF4 import Dataset, date2index, num2date,MFDataset
from mpl_toolkits.basemap import Basemap,addcyclic
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib
from matplotlib.pyplot import *
import os
from scipy import interpolate,signal
import numpy
import numpy as np
numpy.set_printoptions( threshold=500)
import struct
import glob
from matplotlib.patches import Polygon

def xshow(filename, nx, nz):
    f = open(filename, "rb")
    pic = np.zeros((nx, nz))
    for i in range(nx):
        for j in range(nz):
            data = f.read(4)
            elem = struct.unpack("f", data)[0]
            pic[i][j] = elem
    f.close()
    return pic
name=[]

landsat =[]
for nn,ii in enumerate(glob.glob('*.dat')):
    name.append(ii[0:8])  
    landsat.append(xshow(ii,3002,3002))
landsat = np.array(landsat)
landsat1 = ma.masked_array(landsat, np.isnan( landsat) | (landsat<0) | (landsat>1) )

fig, ax = plt.subplots(3,2,figsize=(15,15))
for ii in np.arange(6):
    #这个sca是多图的预定义
  sca(ax[ np.unravel_index(ii,(3,2))])
    m = Basemap(projection='cyl',resolution='c',llcrnrlat=25.2,urcrnrlat=26.1,
                                            llcrnrlon=115.9,urcrnrlon=116.8,lon_0=120.,
                                            llcrnrx=0)

    shp_info = m.readshapefile("CHN_adm_shp\\CHN_adm3",'states',drawbounds=False,linewidth = 0.4) 
    #这里开始画县级底图,只显示范围内的,范围外的填充白色
    for info, shp in zip(m.states_info, m.states):
        proid = info['NAME_3'] 
        if proid == 'Changting':
            poly = Polygon(shp,facecolor='None',edgecolor='b', lw=0.2)
      #通过循环定位了底图所在的位置
      ax[ np.unravel_index(ii,(3,2))].add_patch(poly)
        else:
            poly = Polygon(shp,facecolor='w',edgecolor='w', lw=0.8)
            ax[ np.unravel_index(ii,(3,2))].add_patch(poly)

  #设置每张图的标题
    ax[ np.unravel_index(ii,(3,2))].set_title(name[ii])    
    m.imshow(landsat1[ii,::-1,:],norm=matplotlib.colors.BoundaryNorm( np.linspace(0,1,200), cm.rainbow_r.N),cmap='rainbow_r')

出图效果如下

更进一步,现在想要市级底图和县级底图的叠加。思路和上面类似,通过分别设置,把边界画到一张图上。

代码语言:javascript
复制
from mpl_toolkits.basemap import Basemap
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import maskoceans
from mpl_toolkits.basemap import interp
import numpy.ma as ma
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy
import shapefile

fig, ax = plt.subplots()
rh_deviation=nc.Dataset('rh_201609281230_deviation.nc')
lon=rh_deviation.variables['lon']
lon=np.array(lon)
lat=rh_deviation.variables['lat']
lat=np.array(lat)
rh=rh_deviation.variables['rh']
rh=np.array(rh)
rh_use = np.ma.masked_array(rh, np.isnan( rh) | (rh<-1000)  )
m = Basemap(projection='cyl',resolution='i',llcrnrlon=115.6,llcrnrlat=23,urcrnrlon=120.5,urcrnrlat=28.5,lon_0=120.,lat_0=90)
Lon,Lat = np.meshgrid(lon[:],lat[:])
X,Y = m(Lon,Lat)
#保留福建省境内的县一级边界,福建省边界以外的区域全部填充白色
shp_info3 = m.readshapefile("CHN_adm_shp\\CHN_adm3",'states',drawbounds=False,linewidth = 0.4,zorder=10)
for info, shp in zip(m.states_info, m.states):
    proid = info['NAME_1']  # 可以用notepad打开CHN_adm1.csv文件,可以知道'NAME_1'代表各省的名称
    if proid == 'Fujian':
        poly = Polygon(shp,facecolor='None',edgecolor='b', lw=0.8)
        ax.add_patch(poly)
    else:
        poly = Polygon(shp,facecolor='w',edgecolor='w', lw=0.8)
        ax.add_patch(poly)

#市一级底图
shp_info2 = m.readshapefile("CHN_adm_shp\\CHN_adm2",'states',drawbounds=False,linewidth = 0.4,zorder=10)
for info, shp in zip(m.states_info, m.states):
    proid = info['NAME_1']
    if proid == 'Fujian':
        poly = Polygon(shp,facecolor='None',edgecolor='black', lw=0.8)
        ax.add_patch(poly)
cs=m.contourf(Lon,Lat,rh_use[:,:],np.linspace(-20,20,21),cmap='rainbow')
cbar = plt.colorbar(cs)
plt.show()

出图效果

但是我用的是陆面资料,想把海洋上的值去掉那应该怎么办?那么在上一段代码中45行后加上:

代码语言:javascript
复制
#下面这一段代码块是去掉海洋
vertices = []
codes = []
shp_info = m.readshapefile("CHN_adm_shp\\CHN_adm1",'states',drawbounds=False,linewidth = 0.4,zorder=10)
for info, shp in zip(m.states_info, m.states):
    proid = info['NAME_1']    
    if proid == 'Fujian':        
        for j in np.arange(0, len(shp)):
            vertices.append((shp[j][0], shp[j][1]))
        codes += [Path.MOVETO]
        codes += [Path.LINETO] * (len(shp)-2)
        codes += [Path.CLOSEPOLY]
        clip = Path(vertices, codes)
        clip = PathPatch(clip, transform=ax.transData)
for contour in cs.collections:
    contour.set_clip_path(clip)
##这里结束###############################

效果如下

现在来最终加强版多图的市县边界叠加和海洋掩膜。

代码语言:javascript
复制
from mpl_toolkits.basemap import Basemap
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.pyplot import *
from mpl_toolkits.basemap import maskoceans
from mpl_toolkits.basemap import interp
import numpy.ma as ma
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from mpl_toolkits.axes_grid1 import make_axes_locatable
from osgeo import gdal
import numpy
import shapefile
import struct
import glob
from mpl_toolkits.basemap import Basemap,addcyclic
from netCDF4 import Dataset, date2index, num2date,MFDataset
from matplotlib.patches import Polygon

fig, ax = plt.subplots(3,6,figsize=(15,15))
file=nc.glob('*.nc')
for nn,ii in enumerate(nc.glob('*.nc')):
    rh_deviation=nc.Dataset(ii)
    lon=rh_deviation.variables['lon']
    lon=np.array(lon)
    lat=rh_deviation.variables['lat']
    lat=np.array(lat)
    rh=rh_deviation.variables['rh']
    rh=np.array(rh)
    m = Basemap(projection='cyl',resolution='i',llcrnrlon=115.6,llcrnrlat=23,urcrnrlon=120.5,urcrnrlat=28.5,lon_0=120.,lat_0=90)
    Lon,Lat = np.meshgrid(lon[:],lat[:])
    X,Y = m(Lon,Lat)
    sca(ax[ np.unravel_index(nn,(3,6))]) 
    rh_use = np.ma.masked_array(rh, np.isnan( rh) | (rh<-1000)  )
    cs=m.contourf(Lon,Lat,rh_use,np.linspace(-20,20,21),cmap='rainbow')
#县一级底图    
    shp_info3 = m.readshapefile("CHN_adm_shp\\CHN_adm3",'states',drawbounds=False,linewidth = 0.4,zorder=10)
    for info, shp in zip(m.states_info, m.states):
        proid = info['NAME_1']  # 可以用notepad打开CHN_adm1.csv文件,可以知道'NAME_1'代表各省的名称
        if proid == 'Fujian':
            poly = Polygon(shp,facecolor='None',edgecolor='b', lw=0.3)
            ax[ np.unravel_index(nn,(3,6))].add_patch(poly)
        else:
            poly = Polygon(shp,facecolor='w',edgecolor='w', lw=0.3)
            ax[ np.unravel_index(nn,(3,6))].add_patch(poly)
#市一级底图
    shp_info2 = m.readshapefile("CHN_adm_shp\\CHN_adm2",'states',drawbounds=False,linewidth = 0.4,zorder=10)
    for info, shp in zip(m.states_info, m.states):
        proid = info['NAME_1']
        if proid == 'Fujian':
            poly = Polygon(shp,facecolor='None',edgecolor='black', lw=0.5)
            ax[ np.unravel_index(nn,(3,6))].add_patch(poly)
#去掉海上
    vertices = []
    codes = []    
    shp_info = m.readshapefile("CHN_adm_shp\\CHN_adm1",'states',drawbounds=False,linewidth = 0.4,zorder=10)
    for info, shp in zip(m.states_info, m.states):
        proid = info['NAME_1']      
        if proid == 'Fujian':        
            for j in np.arange(0, len(shp)):
                vertices.append((shp[j][0], shp[j][1]))
            codes += [Path.MOVETO]
            codes += [Path.LINETO] * (len(shp)-2)
            codes += [Path.CLOSEPOLY]
            clip = Path(vertices, codes)
            clip = PathPatch(clip, transform=ax[ np.unravel_index(nn,(3,6))].transData)

    for contour in cs.collections:
        contour.set_clip_path(clip)    
    ax[ np.unravel_index(nn,(3,6))].set_title(name[nn])
fig.colorbar(cs,  ax=ax.ravel().tolist(), orientation='horizontal')
plt.show()

效果如下

另外在前几天MeteoAI中的群对话有这样的信息。

有大佬提问:图中已经利用ETOPO水深数据标记了陆地,但是湖泊还在,只想留下海洋,python有什么合适的方法或者库把湖泊也标记出来

另一位大佬回答,利用from skimage import measure。先把图像二值化,然后找到所有单连通区域,最大单联通区域是海洋通过这种想法把海洋保留陆地去掉。

友情提醒呢,由于本人非地图学专业,很多shp的操作或者提供的代码块也没弄懂,全靠自己摸索。

References

[1]白化:http://bbs.06climate.com/forum.php?mod=viewthread&tid=42437&highlight=%B0%D7%BB%AF

[2]python站点资料插值画图及白化:http://bbs.06climate.com/forum.php?mod=viewthread&tid=57932&highlight=%B0%D7%BB%AF

[3]Clippinga raster with a shapefile:https://basemaptutorial.readthedocs.io/en/latest/clip.html

[4]判断点是不是在shp内:https://blog.csdn.net/liuchengzimozigreat/article/details/93055078

[5]Meshgrid函数的基本用法:https://www.cnblogs.com/lantingg/p/9082333.html

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2019-12-19,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 MeteoAI 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • References
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档