利用SHP裁剪NetCDF数据

网上有很多利用SHP裁剪NC数据的方法，如利用arcgis、arcpy、Qgis、CDO、salem等等，本人也写过一个，在这里再次分享一个比较简单的在Python中利用GDAL裁剪NetCDF数据的方法。

数据准备（关注末尾公众号获得）

一：带有越南行政边界的shapefile

二：netCDF 文件

py库安装

pip install gdal

导入所需要的包

import netCDF4
import numpy as np
from osgeo import gdal,osr,ogr
import matplotlib.pyplot as plt

主函数

#函数：利用shp数据掩膜
#参数：lon&lat经纬度，res像元尺寸，shapefile矢量文件
def makeMask(lon,lat,res,shapefile):
    source_ds = ogr.Open(shapefile)
   
    source_layer = source_ds.GetLayer()
 
    # Create high res raster in memory
    mem_ds = gdal.GetDriverByName('MEM').Create('', lon.size, lat.size, gdal.GDT_Byte)
    mem_ds.SetGeoTransform((lon.min(), res, 0, lat.max(), 0, -res))
    band = mem_ds.GetRasterBand(1)
 
    # Rasterize shapefile to grid
    gdal.RasterizeLayer(mem_ds, [1], source_layer, burn_values=[1])
 
    # Get rasterized shapefile as numpy array
    array = band.ReadAsArray()
 
    # Flush memory file
    mem_ds = None
    band = None
    return array

降水数据展示

# 路径设置
plt.figure(dpi=600)
shapefile = "D:/路径/VNM_adm0.shp"
ncs= "D:/路径/ds_tmys_mmsVN_1976-2005_pr_historical_GFDL-ESM2M.nc"

 
# 读入NC文件
nc = netCDF4.Dataset(ncs,'r')
 
# 得到降水数据
pr = nc.variables['pr'][:]
 
# 展示降水数据
plt.imshow(pr)
plt.show()

 掩膜数据展示

# 得到经纬度信息
lons = nc.variables['lon'][:]

lats = nc.variables['lat'][:]
# 计算像元尺寸
cellsize = lons[:][1] - lons[:][0]

# 创造掩膜数据
mask = makeMask(lons,lats,cellsize,shapefile)
 
# 展示掩膜数据
plt.imshow(mask)
plt.show()

对降水数据进行掩膜

# 对降水数据进行掩膜
precip = np.ma.masked_where(mask==0,pr)
 
plt.imshow(precip)
plt.show()
 
# 输出掩膜后一些计算指标
print (np.min(precip), np.mean(precip), np.max(precip))

 

 这样数据就裁剪好了，个人觉得在处理NC数据时先不要裁剪，在算力允许的情况下整体运算最后裁剪更加便利。

欢迎关注遥感迷公众号，回复矢量裁剪NC获取数据下载链接