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.

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行后加上这一块代码段，于是乡镇边界上去了。

#画出乡镇边界这段为这段代码
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)
#################################################################

出图效果如下

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

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')

出图效果如下

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

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行后加上：

#下面这一段代码块是去掉海洋
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)
##这里结束###############################

效果如下

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

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()

效果如下