修复自制光栅和形状文件叠加不正确的问题

问题描述 投票:0回答:1

我正在为我正在进行的项目的一部分创建巴西南里奥格兰德州的栅格。它们已成功生成,并且一切看起来都很好,除了当我将其绘制在与相同状态的 shapefile 相同的轴上时之外。两者应重叠为相同的大小和相对形状,但 SHP 的绘图较小且上下颠倒。我需要解决这个问题,因为我的下一步是将它们叠加起来进行区域统计。

这是我用来剪辑从中获取数据的初始 netcdf 的代码,然后是我用来从 xarray 数据集中写出栅格的代码,以防在初始创建中出现任何明显的错误。

# Assign bounding lat/lon around RGdS
min_lat, max_lat = -37.0, -23.0
min_lon, max_lon = -61.0, -46.0

# Open dataset, clip
ds = xarray.open_dataset("path")
ds_clipped = ds.sel(lat=slice(min_lat, max_lat), lon=slice(min_lon, max_lon))


# Later on, create a raster from an xarray dataset called GDD_arr
outpath = "path"
with rasterio.open(
    outpath,
    mode="w",
    driver="GTiff",
    height=GDD_arr.shape[1],
    width=GDD_arr.shape[2],
    count=1,
    dtype='float32',
    crs= CRS.from_epsg(4326),
    transform = from_origin(-60.875, -23.125, 0.25, 0.25),
    nodata = -9999) as new_dataset:
  new_dataset.write(GDD_arr[0], 1)

这也是我用来绘制栅格和 SHP 的代码及其外观:

# opening the raster and shapefile
tiff = rasterio.open("GDD_Count_ssp245_2020_tobacco.tif")
shapefile = gpd.read_file("RGdS_AgData.shp")

# set up plot
f, ax = plt.subplots()

# plot DEM
rasterplot.show(np.flipud(tiff.read(1)), ax=ax)

# plot shapefiles
shapefile.plot(ax=ax, facecolor='w', edgecolor='k')
plt.show()

我尝试将栅格的经度转换为 -180 到 180,而不是 0 到 360,但结果仍然看起来像这样,有一个垂直翻转的小 SHP。如果有任何帮助的话,我将在下面添加一张图片,显示它的实际外观,以及栅格和 SHP 的 zip。如果有人知道为什么会发生这种情况或者我可以做些什么来让他们排队,我将非常感激!

数据压缩文件链接

python raster geopandas python-xarray shapefile
1个回答
0
投票

在以下行中:

# plot DEM
rasterplot.show(np.flipud(tiff.read(1)), ax=ax)

您正在将 tiff 像素读取到 numpy 数组,然后简单地在 numpy 数组中绘制值。这样,地理配准信息(变换)就会丢失,因此像素仅被绘制为从 (0,0) 原点开始且像素大小为 1 的笛卡尔数据。

尝试使用rasterio绘图功能,因为它们在绘图时会正确使用地理配准数据。例如:

from rasterio.plot import show

# plot DEM
show(tiff, ax=ax)
© www.soinside.com 2019 - 2024. All rights reserved.