開心市民孫先生

只要是地理專業，就繞不開出圖問題，不論是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地理出图配色及旁门左道的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。