0

各位晚上好。我对 rasterio 中的翘曲选项有疑问。我有两个来自阿尔巴尼亚的栅格数据集。第一个栅格是裁剪的产物，第二个栅格是国家的边界​​。首先，我有以下预测。 在此处输入图像描述

对于第二个栅格，我有以下信息： 在此处输入图像描述

我的目标是使用栅格 2 的空间信息重新投影栅格#1。我尝试使用 rasterio，结果显示以下空间参考：Krassovsky_1942_Transverse_Mercator，并且基准未知。使用以下代码后得到此结果：

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling


#-----------------------------------------------------------------------------  #
#                   OPENING IMAGE WITH THE CORRECT SPATIAL REFERENCE
#-----------------------------------------------------------------------------#
with rasterio.open(r"C:\Users\Roger\Documents\git\ArcpyRecipes\Test_rasterio\test_albania\raster\AL_020m_nat_buffer100m.tif") as proj:
dst_crs = proj.crs
print dst_crs
transform_proj, width_proj, height_proj = calculate

#-----------------------------------------------------------------------------#
#                   OPENING IMAGE WITH TO BE REPROJECTED
#-----------------------------------------------------------------------------#

with rasterio.open(r"C:\Users\Roger\Documents\git\ArcpyRecipes\Test_rasterio\test_albania\results\umprojected\ALB_adm0.tif") as src:
transform, width, height = calculate_default_transform(
    src.crs, dst_crs, src.width, src.height, *src.bounds)
kwargs = src.meta.copy()
kwargs.update({
    "crs": dst_crs,
    "transform": transform,
    "width": width,
    "height": height
})

#---------------------------------------------------------------------------
#                   CREATING THE OUTPUT IMAGE   
#---------------------------------------------------------------------------

with rasterio.open(r"C:\Users\Roger\Documents\git\ArcpyRecipes\Test_rasterio\test_albania\results\projected/Albita_projected_2.tif", "w", **kwargs) as dst:
    for i in range(1, src.count + 1):
        reproject(
            source = rasterio.band(src, i),
            destination=rasterio.band(dst, i),
            src_transform = src.transform,
            src_crs=src.crs,
            dst_transform = transform,
            dst_crs = dst_crs,
            resampling = Resampling.nearest)

我读到可以使用 arcpy.Describe 在 ArcMap 中检索 EPGS 代码，但我想使用模块 rasterio 来完成任务。我认为使用光栅打开图像后的选项（dataset.crs）会给我 CRS 进行投影，但事实并非如此。谁能帮我解决这个问题？

非常感谢

4

0 回答 0