PYTHON地理出图及旁门左道

josephsun

只要是地理专业，就绕不开出图问题，不论是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 第三方库上下载安装。

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

1. import warnings
   
2. warnings.simplefilter('ignore')
   
3. from matplotlib.pylab import rcParams
   
4. rcParams['figure.figsize'] = 14, 14
   
5. import matplotlib.pyplot as plt
   
6. from mpl_toolkits.basemap import Basemap
   
7. import cartopy.crs as ccrs
   
8. import cartopy.mpl.ticker as cticker
   
9. import cartopy.io.shapereader as shpreader
   
10. from osgeo import gdal
    
11. import netCDF4 as nc
    
12. import numpy as np
    
13. import pandas as pd
    
14. import glob

复制代码

读取tif文档

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

复制代码

读取tif函数

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

复制代码

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

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

复制代码

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

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

复制代码

出图样子

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

复制代码

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

1. import geopandas as gp
   
2. import regionmask
   
3. shp = gp.read_file("shp/Uzb.shp")
   
4. mask = regionmask.mask_geopandas(shp, lons, lats)
   
5. def makemask(data):
   
6.     return np.ma.masked_array(data, mask=mask)fig2 = plt.figure()
   
7. proj = ccrs.PlateCarree(central_longitude=lons[int((len(lons)+1)/2)])
   
8. #设置一个圆柱投影坐标，中心经度设置在中间位置
   
9. leftlon, rightlon, lowerlat, upperlat = (55,74.0001,37,46.0001)
   
10. #设置地图边界范围
    
11. lon_formatter = cticker.LongitudeFormatter()
    
12. lat_formatter = cticker.LatitudeFormatter()
    
13. m = Basemap(leftlon,lowerlat,rightlon,upperlat)
    
14. x, y = m(*np.meshgrid(lons, lats))
    
15. f2_ax1 = fig2.add_axes([0.67, 0.8, 1, 1],projection = proj)
    
16. f2_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
    
17. # f2_ax1.set_title('(b)',loc='left',fontsize =15)
    
18. f2_ax1.set_xticks(np.arange(leftlon,rightlon,5), crs=ccrs.PlateCarree())
    
19. f2_ax1.set_yticks(np.arange(lowerlat,upperlat,3), crs=ccrs.PlateCarree())
    
20. f2_ax1.xaxis.set_major_formatter(lon_formatter)
    
21. f2_ax1.yaxis.set_major_formatter(lat_formatter)
    
22. shp = shpreader.Reader("shp/Uzb.shp").geometries()
    
23. f2_ax1.add_geometries(shp, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
    
24. c1 = f2_ax1.contourf(x, y,makemask(readTif(data_list[0])), zorder=0,extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdYlBu)
    
25. plt.tick_params(labelsize=22)
    
26. fig2.savefig('2.png',dpi=300,bbox_inches='tight')
    
27. plt.show()

复制代码

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

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

复制代码

试试出图

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

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

1. plt.rcParams["font.sans-serif"]=["Times New Roman"] #设置字体
   
2. plt.rcParams["axes.unicode_minus"]=False #该语句解决图像中的“-”负号的乱码问题
   
3. fig1 = plt.figure()
   
4. proj = ccrs.PlateCarree(central_longitude=lons[int((len(lons)+1)/2)])
   
5. #设置一个圆柱投影坐标，中心经度设置在中间位置
   
6. leftlon, rightlon, lowerlat, upperlat = (55,74.0001,37,46.0001)
   
7. #设置地图边界范围
   
8. lon_formatter = cticker.LongitudeFormatter()
   
9. lat_formatter = cticker.LatitudeFormatter()
   
10. m = Basemap(leftlon,lowerlat,rightlon,upperlat)
    
11. x, y = m(*np.meshgrid(lons, lats))
    
12. 
    
13. f2_ax1 = fig1.add_axes([0, 1.2, 1, 0.5],projection = proj)
    
14. DrawShp("(a)")
    
15. plt.ylabel('Jan.-Apr. \n', fontsize=45) # 列title
    
16. plt.title('Jan.', fontsize=50) # title
    
17. cbposition=fig1.add_axes([0.9,1.52, 0.01, 0.15])
    
18. DrawTif(data_list[0],0)
    
19. 
    
20. f2_ax1 = fig1.add_axes([1.15, 1.2, 1, 0.5],projection = proj)
    
21. DrawShp("(b)")
    
22. plt.title('Feb.', fontsize=50) # title
    
23. cbposition=fig1.add_axes([2.05,1.52, 0.01, 0.15])
    
24. DrawTif(data_list[1],0)
    
25. 
    
26. f2_ax1 = fig1.add_axes([2.3, 1.2, 1, 0.5],projection = proj)
    
27. DrawShp("(c)")
    
28. plt.title('Mar.', fontsize=50) # title
    
29. cbposition=fig1.add_axes([3.2,1.52, 0.01, 0.15])
    
30. DrawTif(data_list[2],0)
    
31. 
    
32. f2_ax1 = fig1.add_axes([3.45, 1.2, 1, 0.5],projection = proj)
    
33. DrawShp("(d)")
    
34. plt.title('Apr.', fontsize=50) # title
    
35. cbposition=fig1.add_axes([4.35,1.52, 0.01, 0.15])
    
36. DrawTif(data_list[3],0)
    
37. 
    
38. f2_ax1 = fig1.add_axes([0, 0.6, 1, 0.5],projection = proj)
    
39. DrawShp("(e)")
    
40. plt.ylabel('Jan.-Apr. \n', fontsize=45) # 列title
    
41. plt.title('May.', fontsize=50) # title
    
42. cbposition=fig1.add_axes([0.9,0.92, 0.01, 0.15])
    
43. DrawTif(data_list[4],0)
    
44. 
    
45. f2_ax1 = fig1.add_axes([1.15, 0.6, 1, 0.5],projection = proj)
    
46. DrawShp("(f)")
    
47. plt.title('Jun.', fontsize=50) # title
    
48. cbposition=fig1.add_axes([2.05,0.92, 0.01, 0.15])
    
49. DrawTif(data_list[5],0)
    
50. 
    
51. f2_ax1 = fig1.add_axes([2.3, 0.6, 1, 0.5],projection = proj)
    
52. DrawShp("(g)")
    
53. plt.title('Jul.', fontsize=50) # title
    
54. cbposition=fig1.add_axes([3.2,0.92, 0.01, 0.15])
    
55. DrawTif(data_list[6],0)
    
56. 
    
57. f2_ax1 = fig1.add_axes([3.45, 0.6, 1, 0.5],projection = proj)
    
58. DrawShp("(h)")
    
59. plt.title('Aug.', fontsize=50) # title
    
60. cbposition=fig1.add_axes([4.35,0.92, 0.01, 0.15])
    
61. DrawTif(data_list[7],0)
    
62. 
    
63. f2_ax1 = fig1.add_axes([0, 0, 1, 0.5],projection = proj)
    
64. DrawShp("(i)")
    
65. plt.ylabel('Sep.-Oct. \n', fontsize=45) # 列title
    
66. plt.title('Sep.', fontsize=50) # title
    
67. cbposition=fig1.add_axes([0.9,0.32, 0.01, 0.15])
    
68. DrawTif(data_list[8],0)
    
69. 
    
70. f2_ax1 = fig1.add_axes([1.15, 0, 1, 0.5],projection = proj)
    
71. DrawShp("(j)")
    
72. plt.title('Oct.', fontsize=50) # title
    
73. cbposition=fig1.add_axes([2.05,0.32, 0.01, 0.15])
    
74. DrawTif(data_list[9],0)
    
75. 
    
76. f2_ax1 = fig1.add_axes([2.3, 0, 1, 0.5],projection = proj)
    
77. DrawShp("(k)")
    
78. plt.title('Nov.', fontsize=50) # title
    
79. cbposition=fig1.add_axes([3.2,0.32, 0.01, 0.15])
    
80. DrawTif(data_list[10],0)
    
81. 
    
82. f2_ax1 = fig1.add_axes([3.45, 0, 1, 0.5],projection = proj)
    
83. DrawShp("(l)")
    
84. plt.title('Dec.', fontsize=50) # title
    
85. cbposition=fig1.add_axes([4.35,0.32, 0.01, 0.15])
    
86. DrawTif(data_list[11],0)
    
87. 
    
88. fig1.savefig('4.png',dpi=300,bbox_inches='tight')
    
89. plt.show()

复制代码

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

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

josephsun

自己顶哈{:eb502:}{:eb502:}{:eb502:}{:eb502:}{:eb502:}

lzyyy

顶一个

千与千叶

顶一个

流浪的心

顶一个

AllenLu

厉害，顶一个！

WesleyMoob

厉害顶一个！

xiaomaozhao77

非常厉害啊，学习一下

yqcxm

顶一个先！

小妍ivy

非常厉害啊，学习一下