shp文件裁剪多图层nc文件并逐层转为tif的python代码

最新推荐文章于 2024-03-24 11:01:34 发布

PI理论

最新推荐文章于 2024-03-24 11:01:34 发布

阅读量522

点赞数 6

分类专栏： Python 文章标签： python

本文链接：https://blog.csdn.net/weixin_46080995/article/details/135863285

版权

Python 专栏收录该内容

57 篇文章 1 订阅

订阅专栏

本文介绍了如何使用Python库如geopandas和rasterio读取shapefile和NetCDF文件，对黄河流域数据进行空间裁剪，然后按时间序列将结果转换为GeoTIFF格式，以显示在特定坐标系统下特定时间段的辐射数据。

摘要由CSDN通过智能技术生成

import os
import geopandas as gpd
import pandas as pd
import xarray as xr
import rasterio
from matplotlib import pyplot as plt
from rasterio.mask import mask
from rasterio.transform import from_origin
from rasterio.enums import Resampling
from shapely.geometry import mapping

# 尝试用黄河流域的shp文件裁剪nc
# 20240103 09:50 YMJ

# 1. 读取区域形状文件（shp）
shp_path = 'D:\\CMFD\\test\\area.shp'
gdf = gpd.read_file(shp_path)

# 2. 读取NetCDF文件,并定义坐标系
nc_path = 'D:\\CMFD\\test\\lrad_CMFD_V0106_B-01_01mo_010deg_197901-201812.nc'
ds = xr.open_dataset(nc_path)
ds.rio.write_crs("EPSG:4326", inplace=True)

# 3. 根据区域形状文件裁剪NetCDF数据
ds.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=True)
clipped_ds = ds.rio.clip(gdf.geometry.apply(mapping),gdf.crs,drop=True)


# 获取时间信息并将每层保存为TIFF文件
output_directory = 'D:\\CMFD\\test\\output\\'
os.makedirs(output_directory, exist_ok=True)

# 4. 获取时间信息并将每层保存为TIFF文件
for time_index, time_value in enumerate(clipped_ds['time'].values):
    # 提取对应时间层的数据
    data_layer = clipped_ds['lrad'].isel(time=time_index)
    # print(data_layer)

    # 获取时间信息, 保存为TIFF文件
    time_info = pd.to_datetime(ds["time"].isel(time=time_index).values)
    output_path = "D:\\CMFD\\test\\output\\"
    tiff_path = f"{output_path}prec_{time_info.strftime('%Y%m%d')}.tif"
    # print(tiff_path)

    # 获取空间参考和变换信息
    # rasterio.transform.from_origin(west, north, xsize, ysize)
    transform = from_origin(data_layer['lon'].values.min(), data_layer['lat'].values.min(),
                            data_layer['lon'].values[1] - data_layer['lon'].values[0],
                            data_layer['lat'].values[0] - data_layer['lat'].values[1])

    # 写入TIFF文件
    with rasterio.open(tiff_path, 'w', driver='GTiff', height=data_layer.shape[0],
                       width=data_layer.shape[1], count=1, dtype=str(data_layer.dtype),
                       crs='EPSG:4326', transform=transform) as dst:
        dst.write(data_layer.values, 1)

# 关闭NetCDF文件
ds.close()