python 海底地形图 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, 14

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

import cartopy.crs as ccrs
import cartopy.mpl.ticker as cticker
import cartopy.io.shapereader as shpreader

from osgeo import gdal

import netCDF4 as nc

import 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 regionmask

shp = 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_locator
    plt.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) # title

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