PYTHON地理出图配色及旁门左道

开心市民孙先生

只要是地理专业，就绕不开出图问题，不论是arcgis、supermap的基本配色都比较固定，容易造成审美疲劳，Qgis自定义是个加分项，但基本功能还是不足以满足空间分析。关于批量显示tif问题，不论是哪个软件都较为繁琐，需要单独添加，花费大量时间。

本章将基于Python批量化出图问题，涉及一些地理库及地理掩膜功能的操作介绍，关于部分库使用及版本问题，可以公众号留言交流。大多数的气象数据都是NC数据文件，本章不涉及NC转TIF过程，后续会专门出一章NC批量转TIF并计算年降水保存为tif，本章节将从读取tif作为示例，有关批量显示制图技巧的详情参见原文链接。

涉及到的出图库包括GDAL、matplotlib、cartopy、Basemap等，对于库的介绍请参考官方文档，在此不多赘述。如果pip装不上可以从 https://www.lfd.uci.edu/~gohlke/pythonlibs/#pylibtiff 第三方库上下载安装。

首先导入基本应用库[不一定都会用到]

import warnings
warnings.simplefilter('ignore')from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 14, 14import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemapimport cartopy.crs as ccrs
import cartopy.mpl.ticker as cticker
import cartopy.io.shapereader as shpreaderfrom osgeo import gdalimport netCDF4 as ncimport numpy as np
import pandas as pd
import glob

读取tif文档

# 读取所有nc数据
Input_folder = './data'
data_list = glob.glob(Input_folder + '/*.tif')
data_list

读取tif函数

def readTif(fileName):dataset = gdal.Open(fileName)im_width = dataset.RasterXSize #栅格矩阵的列数im_height = dataset.RasterYSize #栅格矩阵的行数print(im_width, im_height)im_data = dataset.ReadAsArray(0,0,im_width,im_height)#获取数据im =  im_data[0:im_height,0:im_width]#获取蓝波段return im

tif基础信息[批量出图需要有一个参考]

dataset = gdal.Open(data_list[0])
im_width = dataset.RasterXSize #栅格矩阵的列数
im_height = dataset.RasterYSize #栅格矩阵的行数im_geotrans = dataset.GetGeoTransform()#获取仿射矩阵信息
im_proj = dataset.GetProjection()#获取投影信息
# im =  im_data[0:im_height,0:im_width]#获取蓝波段
im_geotrans

手动设置经纬度[其实有更好方法，读取文件基本属性，本文选用numpy造数组形式]

lons = np.arange(im_geotrans[0], im_geotrans[0]+im_geotrans[1]*im_width, im_geotrans[1])
lats = np.arange(im_geotrans[3], im_geotrans[3]+im_geotrans[5]*im_height, im_geotrans[5])

出图样子

fig2 = plt.figure()
proj = ccrs.PlateCarree(central_longitude=lons[int((len(lons)+1)/2)])
#设置一个圆柱投影坐标，中心经度设置在中间位置
leftlon, rightlon, lowerlat, upperlat = (55,74.0001,37,46.0001)
#设置地图边界范围
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
m = Basemap(leftlon,lowerlat,rightlon,upperlat)
x, y = m(*np.meshgrid(lons, lats))f2_ax1 = fig2.add_axes([0.67, 0.8, 1, 1],projection = proj)
f2_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
# f2_ax1.set_title('(b)',loc='left',fontsize =15)
f2_ax1.set_xticks(np.arange(leftlon,rightlon,5), crs=ccrs.PlateCarree())
f2_ax1.set_yticks(np.arange(lowerlat,upperlat,3), crs=ccrs.PlateCarree())
f2_ax1.xaxis.set_major_formatter(lon_formatter)
f2_ax1.yaxis.set_major_formatter(lat_formatter)shp = shpreader.Reader("shp/Uzb.shp").geometries()
f2_ax1.add_geometries(shp, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)c1 = f2_ax1.contourf(x, y, readTif(data_list[0]), zorder=0,extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdYlBu)
plt.tick_params(labelsize=22)
fig2.savefig('scatter1.png',dpi=300,bbox_inches='tight')
plt.show()

看起来还不错，但是如果批量出图就需要考虑每张图片的位置和图例问题，本文已最极端情况，每张图片图例都需要重新出，而且仅显示shp内数据，在空白地方显示图例，进行说明。

import geopandas as gp
import regionmaskshp = gp.read_file("shp/Uzb.shp")
mask = regionmask.mask_geopandas(shp, lons, lats)def makemask(data):return np.ma.masked_array(data, mask=mask)

fig2 = plt.figure()
proj = ccrs.PlateCarree(central_longitude=lons[int((len(lons)+1)/2)])
#设置一个圆柱投影坐标，中心经度设置在中间位置
leftlon, rightlon, lowerlat, upperlat = (55,74.0001,37,46.0001)
#设置地图边界范围
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
m = Basemap(leftlon,lowerlat,rightlon,upperlat)
x, y = m(*np.meshgrid(lons, lats))f2_ax1 = fig2.add_axes([0.67, 0.8, 1, 1],projection = proj)
f2_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
# f2_ax1.set_title('(b)',loc='left',fontsize =15)
f2_ax1.set_xticks(np.arange(leftlon,rightlon,5), crs=ccrs.PlateCarree())
f2_ax1.set_yticks(np.arange(lowerlat,upperlat,3), crs=ccrs.PlateCarree())
f2_ax1.xaxis.set_major_formatter(lon_formatter)
f2_ax1.yaxis.set_major_formatter(lat_formatter)shp = shpreader.Reader("shp/Uzb.shp").geometries()
f2_ax1.add_geometries(shp, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)c1 = f2_ax1.contourf(x, y,makemask(readTif(data_list[0])), zorder=0,extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdYlBu)
plt.tick_params(labelsize=22)
fig2.savefig('2.png',dpi=300,bbox_inches='tight')
plt.show()

为了批量显示图，减少代码量，需要建立出图函数。

# 画坐标系和shp
def DrawShp(label, data):f2_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())f2_ax1.set_title(label,loc='left',fontsize =50)f2_ax1.set_xticks(np.arange(leftlon,rightlon,5), crs=ccrs.PlateCarree())f2_ax1.set_yticks(np.arange(lowerlat,upperlat,3), crs=ccrs.PlateCarree())f2_ax1.xaxis.set_major_formatter(lon_formatter)f2_ax1.yaxis.set_major_formatter(lat_formatter)shp = shpreader.Reader("shp/Uzb.shp").geometries()f2_ax1.add_geometries(shp, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)plt.tick_params(labelsize=40)# 画tif
def DrawTif(data, m):clevs = np.linspace(np.nanmin(makemask(readTif(data_list[0]))),np.nanmax(makemask(readTif(data_list[0]))),11)c1 = f2_ax1.contourf(x, y,makemask(readTif(data)),clevs, zorder=0,extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdYlBu)tick_locator = ticker.MaxNLocator(nbins=2)if m == 1:cb= fig1.colorbar(c1,cax=cbposition, format='%d')# colorbar上的刻度值个数else: cb= fig1.colorbar(c1,cax=cbposition, format='%.2f')# colorbar上的刻度值个数cb.locator = tick_locatorplt.tick_params(labelsize=35)cb.set_ticks([np.nanmin(makemask(readTif(data))),np.nanmax(makemask(readTif(data)))])cb.outline.set_visible(False)cb.update_ticks()

试试出图

plt.rcParams["font.sans-serif"]=["Times New Roman"] #设置字体fig1 = plt.figure()
proj = ccrs.PlateCarree(central_longitude=lons[int((len(lons)+1)/2)])
#设置一个圆柱投影坐标，中心经度设置在中间位置
leftlon, rightlon, lowerlat, upperlat = (55,74.0001,37,46.0001)
#设置地图边界范围
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
m = Basemap(leftlon,lowerlat,rightlon,upperlat)
x, y = m(*np.meshgrid(lons, lats))# 设置图片位置[横轴坐标，纵轴坐标，长，宽]
# 如果想要最好的长宽比，可以lens(lons)/lens(lats)查看比例
f2_ax1 = fig1.add_axes([0, 1.2, 1, 0.5],projection = proj)
DrawShp("(a)", data_list[0])
plt.ylabel('Jan.-Apr. \n', fontsize=45) # 列title
plt.title('Jan.', fontsize=50) # titlecbposition=fig1.add_axes([0.9,1.52, 0.01, 0.15])
DrawTif(data_list[0],0)
fig2.savefig('3.png',dpi=300,bbox_inches='tight')
plt.show()

关于批量出图，重点在于考虑两个数字

第一图片位置[0, 1.2, 1, 0.5]

第一坐标轴位置[0.9,1.52, 0.01, 0.15]

第二图片位置[1.15, 1.2, 1, 0.5]

第二坐标轴位置[2.05,1.52, 0.01, 0.15]

按照规律即可算出每张图片之间间距，这对于批量出图需要不断尝试间距，寻找最合适间距。

plt.rcParams["font.sans-serif"]=["Times New Roman"] #设置字体
plt.rcParams["axes.unicode_minus"]=False #该语句解决图像中的“-”负号的乱码问题fig1 = plt.figure()
proj = ccrs.PlateCarree(central_longitude=lons[int((len(lons)+1)/2)])
#设置一个圆柱投影坐标，中心经度设置在中间位置
leftlon, rightlon, lowerlat, upperlat = (55,74.0001,37,46.0001)
#设置地图边界范围
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
m = Basemap(leftlon,lowerlat,rightlon,upperlat)
x, y = m(*np.meshgrid(lons, lats))f2_ax1 = fig1.add_axes([0, 1.2, 1, 0.5],projection = proj)
DrawShp("(a)")
plt.ylabel('Jan.-Apr. \n', fontsize=45) # 列title
plt.title('Jan.', fontsize=50) # title
cbposition=fig1.add_axes([0.9,1.52, 0.01, 0.15])
DrawTif(data_list[0],0)f2_ax1 = fig1.add_axes([1.15, 1.2, 1, 0.5],projection = proj)
DrawShp("(b)")
plt.title('Feb.', fontsize=50) # title
cbposition=fig1.add_axes([2.05,1.52, 0.01, 0.15])
DrawTif(data_list[1],0)f2_ax1 = fig1.add_axes([2.3, 1.2, 1, 0.5],projection = proj)
DrawShp("(c)")
plt.title('Mar.', fontsize=50) # title
cbposition=fig1.add_axes([3.2,1.52, 0.01, 0.15])
DrawTif(data_list[2],0)f2_ax1 = fig1.add_axes([3.45, 1.2, 1, 0.5],projection = proj)
DrawShp("(d)")
plt.title('Apr.', fontsize=50) # title
cbposition=fig1.add_axes([4.35,1.52, 0.01, 0.15])
DrawTif(data_list[3],0)f2_ax1 = fig1.add_axes([0, 0.6, 1, 0.5],projection = proj)
DrawShp("(e)")
plt.ylabel('Jan.-Apr. \n', fontsize=45) # 列title
plt.title('May.', fontsize=50) # title
cbposition=fig1.add_axes([0.9,0.92, 0.01, 0.15])
DrawTif(data_list[4],0)f2_ax1 = fig1.add_axes([1.15, 0.6, 1, 0.5],projection = proj)
DrawShp("(f)")
plt.title('Jun.', fontsize=50) # title
cbposition=fig1.add_axes([2.05,0.92, 0.01, 0.15])
DrawTif(data_list[5],0)f2_ax1 = fig1.add_axes([2.3, 0.6, 1, 0.5],projection = proj)
DrawShp("(g)")
plt.title('Jul.', fontsize=50) # title
cbposition=fig1.add_axes([3.2,0.92, 0.01, 0.15])
DrawTif(data_list[6],0)f2_ax1 = fig1.add_axes([3.45, 0.6, 1, 0.5],projection = proj)
DrawShp("(h)")
plt.title('Aug.', fontsize=50) # title
cbposition=fig1.add_axes([4.35,0.92, 0.01, 0.15])
DrawTif(data_list[7],0)f2_ax1 = fig1.add_axes([0, 0, 1, 0.5],projection = proj)
DrawShp("(i)")
plt.ylabel('Sep.-Oct. \n', fontsize=45) # 列title
plt.title('Sep.', fontsize=50) # title
cbposition=fig1.add_axes([0.9,0.32, 0.01, 0.15])
DrawTif(data_list[8],0)f2_ax1 = fig1.add_axes([1.15, 0, 1, 0.5],projection = proj)
DrawShp("(j)")
plt.title('Oct.', fontsize=50) # title
cbposition=fig1.add_axes([2.05,0.32, 0.01, 0.15])
DrawTif(data_list[9],0)f2_ax1 = fig1.add_axes([2.3, 0, 1, 0.5],projection = proj)
DrawShp("(k)")
plt.title('Nov.', fontsize=50) # title
cbposition=fig1.add_axes([3.2,0.32, 0.01, 0.15])
DrawTif(data_list[10],0)f2_ax1 = fig1.add_axes([3.45, 0, 1, 0.5],projection = proj)
DrawShp("(l)")
plt.title('Dec.', fontsize=50) # title
cbposition=fig1.add_axes([4.35,0.32, 0.01, 0.15])
DrawTif(data_list[11],0)fig1.savefig('4.png',dpi=300,bbox_inches='tight')
plt.show()

具体实验操作已上传github上，可以下载尝试，麻烦给个星星哈！~

https://github.com/JosephSun-6/Climate-Variability

PYTHON地理出图配色及旁门左道相关推荐

Python地理做图——学习笔记
Python地理做图--学习笔记 GMT 绘制海岸线绘制地形并叠加海岸线地理信息数据格式在线转换网址适用OSGEO4w可以实现tif转nc,转grd 绘制grd和nc 除了投影方式-X, gmt ...
arcgis10.7利用python批量出图
需求背景现有张一张二张三等100余户宗地shp文件,需要以谷歌影像为底图绘制各户宗地范围图文件准备包含各户shp的gdb文件制作出图模板制作出图模板出图模板设置成果图的标题元素名称为ma ...
python找出图中所有闭合环_求图中的所有闭合环
NetworkX是一个流行的Python包,用于处理许多科学Python发行版中包含的图形.它包括一些计算图圈的算法.尤其是,^{}会回答你的问题.在这种方法的一个警告是必须将图转换为有向图.这意味 ...
python找出图中所有闭合环_这可能是史上最全的 Python 算法集（建议收藏）
△蓝字可关注并标星 -数据分析展示就用DataHunter- 导读:本文是一些机器人算法(特别是自动导航算法)的Python代码合集.其主要特点有以下三点: 选择了在实践中广泛应用的算法: 依赖最少: ...
python箱线图配色_python绘制箱线图
三种方式绘制箱线图 #第一种:直接使用自带的箱线图函数 import pandas as pd import matplotlib.pyplot as plt data=pd.read_excel(& ...
python批量出图
import glob import numpy as np import matplotlib.pyplot as plt import cartopy.io.shapereader as shpr ...
全国多地新冠病例0增长，教你用Python画出当下疫情最火玫瑰图！
CDA数据分析师出品近日,新冠肺炎防控成果的好消息不断. 今天我们聊聊,惊艳的疫情直观图. 据国家卫健委数据统计, 截止至3月10日24时,31省区市累计治愈出院病历超6万,达到61475例. 3 ...
python画人民大会堂_太震撼了，我用python画出全北京的公交线路动图
原标题:太震撼了,我用python画出全北京的公交线路动图今天教大家用pyecharts制作北京市公交线路动态图,这应该是全网唯一一篇能正常运行的教程一.获取百度秘钥首先,本项目需要引用百度地图 ...
arcgis批量出图python代码_【GIS进阶】ArcGIS批量出图_定义出图
今天的文章是浩哥投稿!!! 下图是我欢呼雀跃的样子~~~~~ 本文亮点: 所有步骤都是用ArcGIS中各种工具和软件操作组合,未使用Arcpy与Python等需要使用代码的工具! 这次的这个批量出图又 ...

PYTHON地理出图配色及旁门左道

PYTHON地理出图配色及旁门左道相关推荐

最新文章

热门文章