# Python中关于底图的操作

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

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

Talk is cheap, show me your code.

```from numpy import *
from netCDF4 import Dataset, date2index, num2date,MFDataset
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):
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后，剩下的区域全部填成白色
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)
else:
poly = Polygon(shp,facecolor='w',edgecolor='w', lw=0.8)
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()```

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

#这里开始画县级底图，只显示范围内的，范围外的填充白色
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)
#通过循环定位了底图所在的位置
else:
poly = Polygon(shp,facecolor='w',edgecolor='w', lw=0.8)

#设置每张图的标题
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 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)
#保留福建省境内的县一级边界，福建省边界以外的区域全部填充白色
for info, shp in zip(m.states_info, m.states):
if proid == 'Fujian':
poly = Polygon(shp,facecolor='None',edgecolor='b', lw=0.8)
else:
poly = Polygon(shp,facecolor='w',edgecolor='w', lw=0.8)

#市一级底图
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)
cs=m.contourf(Lon,Lat,rh_use[:,:],np.linspace(-20,20,21),cmap='rainbow')
cbar = plt.colorbar(cs)
plt.show()```

```#下面这一段代码块是去掉海洋
vertices = []
codes = []
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 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 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')
#县一级底图
for info, shp in zip(m.states_info, m.states):
if proid == 'Fujian':
poly = Polygon(shp,facecolor='None',edgecolor='b', lw=0.3)
else:
poly = Polygon(shp,facecolor='w',edgecolor='w', lw=0.3)
#市一级底图
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)
#去掉海上
vertices = []
codes = []
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()```

### References

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

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

0 条评论

• ### Metpy新版功能下载TLnP图设置

用于天气绘图的Metpy包更新（0.8版本）了，他们要逐渐抛弃Python2.X，转到Python>=3.6的版本上。所以，之前（越2018年6月以前，0.7版...

• ### Keras系列 （4）LSTM的返回序列和返回状态的区别

长期短期记忆（LSTM）是由三个内部闸(internal gates)所构建成的循环神经网络(recurrent neuralnetwork)。

• ### Numpy基本用法介绍

我们在以前的文章中已经介绍了如何安装python及其python的一些特性，现在将介绍数据分析过程中经常用到的Numpy库。

• ### Python炫技操作：花式导包的八种方法

__import__ 函数可用于导入模块，import 语句也会调用函数。其定义为：

• ### spring-boot 速成(11) - 单元测试

一、添加依赖项： testCompile 'org.springframework.boot:spring-boot-starter-test:1.5.2.RE...

• ### Angular 使用 RxJS 优化处理多个Http请求

注意：上面的this.http.get... 处理HTTP最好放到单独的Service文件中，再注入到Component。这里为了演示没有这么做。

• ### Android之网络摄像头

实现的功能就是两个手机在一个局域网内可以互相观看对方的摄像头图像,当然如果都是连接公网那么就能远程互看了,,,,和视频聊天差不多,,不过没有声音,,,,,,,,...

• ### JPA关联关系表中加其他字段

JPA是Java Persistence API的简称，中文名Java持久层API，是JDK 5.0注解或XML描述对象－关系表的映射关系，并将运行期的实体[对...

• ### Python数据处理从零开始----第四章（可视化）①③多变量绘图目录

科研工作中我们经常需要把每两个变量之间的关系计算， 然后可以得到一个相关矩阵。 如果两个变量的变化趋势一样， 那么这个值就会大于零， 表明连个变量正相关，值越大...

• ### 集成 SpringBoot 2.3.2 + Shiro 1.5.3 + jwt (无状态)

shiro 集成 jwt 需要禁用 session, 服务器就不用维护用户的状态, 做到无状态调用

MeteoAI