您现在的位置是：首页 >技术交流 >python批量裁剪tif nc并可视化网站首页技术交流

引言

对于空间数据，我们感兴趣的往往是其中的某一部分，对于不需要的部分需要做一些掩膜(Mask)。 比如降水、气温这种数据往往是粗分辨率全球的，如CRU，而且他的存储方式是nc。
但是往往我们只需要某个区域的数据：
把大量的全球nc转为需要的tif时，会造成空间浪费和转换缓慢，这时候就需要先进行裁剪，而且要裁剪nc

Python剪裁nc文件

首先介绍 批量裁剪nc文件

本文使用新一代xarray方法，相比于netCDF4更方便、简洁和高效

salem解决mask

Mask的方法其实很多，xarray的salem只需要一个shp文件就可以搞定所有问题。 salem是xarray的扩展包，集成了一些地球科学数据处理的小工具，其中.roi函数可以根据shp文件提取感兴趣的区域。

```bash
# -*- coding: utf-8 -*-
"""
Created on Sat Apr 15 17:24:53 2023

@author: CH
"""
#%%   原文https://mp.weixin.qq.com/s/AJpR7R7veTXvS7yALl5Mgw
#%%读取边界 并查看

# filepath = os.path.join('data', 'D:/C/笔记5——基于cartopy绘制contour并对中国地区进行白化（包含南海）/china_shp')
# data = gpd.read_file(filepath)
import xarray
import geopandas as gpd
from geopandas import GeoDataFrame, read_file
import matplotlib.pyplot as plt
import salem
shp_dir='D:/data/shp/最新2021年全国行政区划/省.shp'
shpfile=gpd.read_file(shp_dir)
fig, ax = plt.subplots(figsize=(12, 12))
ax = shpfile.plot(ax=ax, color='w', edgecolor='r', lw=1)

'''

#%%#%%读取nc
file_name='D:/data/xiangguan_data/air.mon.mean.nc'
ds=xarray.open_dataset(file_name)
print(ds)

```cpp

# 删除 level 维度
# ds = ds.drop('level')

ds = ds.isel(level=0)
print(ds)

#%%
ds['air'].isel(time=0).salem.quick_map()#0时刻快速可视化

```dart
#%%
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
# 导入Cartopy专门提供的经纬度的Formatter
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.ticker as mticker
# 导入底图包
import cartopy.feature as cfeature

import cartopy.io.shapereader as shpreader

data = ds['air']#三个维度
data = np.flipud(data[0, :, :])#time变没了


def add_shp(ax, **kwargs):

    proj = ccrs.PlateCarree()
    reader = shpreader.Reader('D:/data/shp/china_shp/country1.shp')
    provinces = reader.geometries()
    ax.add_geometries(provinces, proj, **kwargs)
    reader.close()
    # reader.plot   
fig = plt.figure()
proj = ccrs.PlateCarree()#ccrs.Robinson()
ax = fig.add_subplot(111, projection=proj)

# 设置经纬度范围,限定为中国
# 注意指定crs关键字,否则范围不一定完全准确
extents = [-180, 180, -90, 90]
ax.set_extent(extents, crs=proj)
# 添加各种特征
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.LAND, edgecolor='black')
#ax.add_feature(cfeature.LAKES, edgecolor='black')
#ax.add_feature(cfeature.RIVERS)
#ax.add_feature(cfeature.BORDERS)
# 添加网格线
ax.gridlines(linestyle='--')

# 设置大刻度和小刻度
ax.set_xticks(np.arange(-180, 180 + 60, 60), crs=proj)
ax.set_xticks(np.arange(-180, 180 + 30, 30), minor=True, crs=proj)
ax.set_yticks(np.arange(-90, 90 + 30, 30), crs=proj)
ax.set_yticks(np.arange(-90, 90 + 15, 15), minor=True, crs=proj)

# 利用Formatter格式化刻度标签
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())

# 添加底图
long = np.linspace(-180, 180, 144); lat = np.linspace(-90, 90, 73)

im = ax.contourf(
    long, lat[::-1], data,
     cmap='RdYlBu_r',
    extend='both', alpha=0.8
)# 
cbar = fig.colorbar(
    im, ax=ax, shrink=0.9, pad=0.1, orientation='horizontal'
)

# 设置显示范围
cbar.set_clim(vmin=-40, vmax=40)

cbar.ax.tick_params(labelsize='small')
ax.set_title('Annul air (mm)')

# 添加矢量
add_shp(ax, lw=1, ec='k', fc='none')
    
plt.show()

只需要用ds.salem.roi(shape=shpfile)即可进行裁剪

#%%裁剪 
maskregion=ds.salem.roi(shape=shpfile)
#
temp = maskregion['air'].isel(time=0)
 temp.salem.quick_map()
#%裁剪后的结果即为maskregion，可以转tif，或者写为nc
maskregion.to_netcdf(path='D:/data/chiana_air.nc')

#%%最后将裁剪结果可视化：
data = maskregion['air']
data = np.flipud(data[0, :, :])
data1 = data[14:30, 100:128]
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
# 导入Cartopy专门提供的经纬度的Formatter
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.ticker as mticker
# 导入底图包
import cartopy.feature as cfeature

import cartopy.io.shapereader as shpreader
def add_shp(ax, **kwargs):
    '''
    在地图上画出中国省界的shapefile.

    Parameters
    ----------
    ax : GeoAxes
        目标地图.

    **kwargs
        绘制shape时用到的参数.例如linewidth,edgecolor和facecolor等.
    '''
    proj = ccrs.PlateCarree()
    reader = shpreader.Reader('D:/data/shp/最新2021年全国行政区划/省.shp')
    provinces = reader.geometries()
    ax.add_geometries(provinces, proj, **kwargs)
    reader.close()
    
fig = plt.figure()
proj = ccrs.PlateCarree()
ax = fig.add_subplot(111, projection=proj)

# 设置经纬度范围,限定为中国
# 注意指定crs关键字,否则范围不一定完全准确
extents = [70, 140, 15, 55]
ax.set_extent(extents, crs=proj)
# 添加各种特征
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.LAND, edgecolor='black')
#ax.add_feature(cfeature.LAKES, edgecolor='black')
#ax.add_feature(cfeature.RIVERS)
#ax.add_feature(cfeature.BORDERS)
# 添加网格线
ax.gridlines(linestyle='--')

# 设置大刻度和小刻度
ax.set_xticks(np.arange(70, 140 + 20, 20), crs=proj)
ax.set_xticks(np.arange(70, 140 + 10, 10), minor=True, crs=proj)
ax.set_yticks(np.arange(15, 55 + 20, 20), crs=proj)
ax.set_yticks(np.arange(15, 55 + 10, 10), minor=True, crs=proj)

# 利用Formatter格式化刻度标签
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())

# 添加底图
long = np.linspace(70, 140, 28); lat = np.linspace(15, 55, 16)

im = ax.contourf(
    long, lat[::-1], data1,levels=np.linspace(-40, 40, 10),
    cmap='RdYlBu_r',
    extend='both', alpha=0.8,vmin=-40, vmax=40
)
cbar = fig.colorbar(
    im, ax=ax, shrink=0.9, pad=0.1, orientation='horizontal'
)

# 设置显示范围
# ScalarMappable.set_clim(vmin=-40, vmax=40)
cbar.ax.tick_params(labelsize='small')
ax.set_title('Annul air ')

# 添加矢量
add_shp(ax, lw=1, ec='k', fc='none')
    
plt.show()

#%%裁剪tiff
from osgeo import gdal,osr,ogr
import os
import glob

def maskTiff(data,Output_folder):
    
    out_tif_name = Output_folder + '\'+ data.split('\')[-1].split('.')[0] + '.tif'
    ds = gdal.Warp(out_tif_name,
                  data,
                  format = 'GTiff',
                  cutlineDSName = input_shape,      
                  cutlineWhere="FIELD = 'whatever'",
                  outputBounds = [70, 15, 140, 55],
                  dstNodata = 0)

Input_folder = 'D:/data/free_study/tif'
Output_folder = 'D:/data/free_study/mask'
input_shape = 'D:/data/shp/最新2021年全国行政区划/省.shp'

data_list = glob.glob(Input_folder + '\*.tif')
for i in range(len(data_list)):
    data = data_list[i]
    maskTiff(data, Output_folder)
    print(str(i) + ': ' + data)

print('ok')

‘’’
这段代码用于对给定的一组 TIFF 格式影像进行裁剪、掩膜操作，并将结果保存为新的 TIFF 格式文件。

具体来说，该代码利用了 osgeo 库中的 gdal 模块，通过 gdal.Warp() 函数实现了对单个 TIFF 影像的裁剪和掩膜操作。其中，data 和 Output_folder 分别表示输入的 TIFF 影像路径和输出掩膜后的 TIFF 影像的保存路径；out_tif_name 表示输出掩膜后的 TIFF 影像的名称，是根据 Output_folder 和输入文件名称生成的；cutlineDSName 和 cutlineWhere 参数分别指定了掩膜形状文件和掩膜条件，用于在裁剪过程中限制输出范围；outputBounds 则指定了输出的空间范围，以地理坐标系值 [xmin, ymin, xmax, ymax] 的方式指定（此代码中的范围对应中国大陆地区）；dstNodata 参数则指定了输出图像中的 nodata 值。

为了对整个文件夹中的 TIFF 影像进行批量处理，该代码利用了 glob() 函数获取输入影像路径，并通过 for 循环遍历每个影像并调用 maskTiff() 函数进行裁剪和掩膜操作。在进行文件处理时，程序会打印出当前处理的文件序号和文件路径，以方便用户进行监控。最终输出结果为裁剪和掩膜后的 TIFF 影像，保存在指定的输出文件夹中。