如何在 python 中使用 shapefile 重叠（mach）和裁剪雷达图像

Question

我正在徘徊是否有人可以帮助我解决以下问题。在此之前，我真的很感激。我有一个 shapefile 和一个雷达图像。但是它们在不同的坐标系中。雷达图像是一个覆盖部分 shapefile 的 NetCDF 文件。我想知道如何将它们重叠（匹配）在一起并从 shapefile 中剪切雷达图像。 shapefile和雷达图像的坐标系如下图

shapefile 的坐标系：离开 2126960.000000 米，右 2938920.000000 米，顶部 2825210.000000 米，底部 2262400.000000 米，投影坐标系：GDA_1994_Lambert_Conformal_Conic，投影：Lambert_Conformal_Conic， False_Easting: 2500000.00000000, False_Northing: 2500000.00000000, Central_Meridian: 145.00000000, Standard_Parallel_1：-36.00000000， Standard_Parallel_2：-38.00000000， Latitude_Of_Origin：-37.00000000，线性单位：米，地理坐标系：GCS_GDA_1994，基准：D_GDA_1994，本初子午线：格林威治， Angular单位：度，

雷达图像的坐标系，y,y_bound,x,x_bound, , int8项目（）， grid_mapping_name: albers_conical_equal_area, standard_parallel: [-36.3 -39.4], longitude_of_central_meridian: 144.752, latitude_of_projection_origin：-37.852， false_easting: 0.0, false_northing: 0.0, semi_major_axis: 6378137.0, semi_minor_axis: 6356752.31414, 无限维度：当前形状 = (), 加油，

是 , float64 y(y), standard_name: projection_y_coordinate, 单位：公里，界限：y_bounds，无限维度：当前形状 = (512,),

Y_bound, float64 y_bounds(y, n2), 无限维度：当前形状 = (512, 2), 填满，加油，

X, float64 x(x), standard_name: projection_x_coordinate, 单位：公里，界限：x_bounds，无限维度：当前形状 = (512,), 加油，

X_bound, float64 x_bounds(x, n2), 无限维度：当前形状 = (512, 2), 填满

Answer 1

如果您需要帮助将不同坐标系的几何图形转换为通用坐标系，您可以使用 pyproj 等。您首先需要找到定义原始数据的坐标系的 epsg 代码。您可以在此处找到它们 https://spatialreference.org/ref/?search=lambert+conformal+conic&srtext=Search。在下面的示例中，我已经尝试使用正确的，但请仔细检查。

from pyproj import CRS
from shapely.geometry import box
from shapely.ops import transform


crs_shapefile = CRS("EPSG:3112")
crs_radar = CRS("EPSG:3577")

transformer = pyproj.Transformer.from_crs(
            crs_from=crs_shapefile,
            crs_to=crs_radar,
        )

example_geometry = box(0,0,1,1)
transformed_box = transform(func=transformer.transform, geom=example_geometry)

如何在 python 中使用 shapefile 重叠（mach）和裁剪雷达图像

How to overlap (mach) and clip a radar image with a shapefile in python

python

numpy

gdal

netcdf

shapely