专栏首页MeteoAIPython中关于底图的操作

Python中关于底图的操作

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

前言

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

效果如下

另外在前几天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

本文分享自微信公众号 - MeteoAI(meteoai),作者:书呆的边缘,小陈

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-12-19

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

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

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

    zhangqibot
  • Keras系列 (4)LSTM的返回序列和返回状态的区别

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

    zhangqibot
  • Numpy基本用法介绍

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

    zhangqibot
  • Python炫技操作:花式导包的八种方法

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

    Python进阶者
  • spring-boot 速成(11) - 单元测试

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

    菩提树下的杨过
  • Angular 使用 RxJS 优化处理多个Http请求

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

    mafeifan
  • Android之网络摄像头

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

    杨奉武
  • JPA关联关系表中加其他字段

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

    用户5166330
  • Python数据处理从零开始----第四章(可视化)①③多变量绘图目录

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

    用户1359560
  • 集成 SpringBoot 2.3.2 + Shiro 1.5.3 + jwt (无状态)

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

    北漂的我

扫码关注云+社区

领取腾讯云代金券