使用不同的投影将经纬度转换为行列时出错

Question

我遵循了这个很好的答案，但是在使用带有 EPSG:21781+5728 的 dem 时我得到了不正确的结果。我正在使用的 dem 可以在here.

下载

我的代码：

# latlong to rowcol function
def latlong_to_rowcol(dem_transform, point_y, point_x):
    topleft_y = dem_transform[3]
    topleft_x = dem_transform[0]
    y_distance = dem_transform[5]
    x_distance = dem_transform[1]
    row = -round((topleft_y - point_y) / y_distance)
    col = round((point_x - topleft_x) / x_distance)
    return row,col

source = osr.SpatialReference()
source.ImportFromEPSG(4326)
target = osr.SpatialReference(wkt=dem.GetProjection())
transform = osr.CoordinateTransformation(source, target)

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(long, lat)
point.Transform(transform)

row,col = latlong_to_rowcol(dem.GetGeoTransform(), point.GetY(), point.GetX())

print("latlong\t\t", lat, long) # 46.911999 7.4520001
print("after transform\t", point.GetY(), point.GetX())
print("rowcol\t\t", row,col)
print("dem.shape\t", demarray.shape)

印花：

latlong      46.911999 7.4520001
after transform  -3040610.7184233973 5670000.516904108
rowcol       133705 207601
dem.shape    (9121, 15401)

如您所见，生成的行和列超出了 dem 形状的范围。

出于好奇，我将 dem 转换为 EPSG:4326 并在新的 dem 上运行代码。上面的代码有效！！！为什么这适用于 EPSG:4326 dem 但不适用于 EPSG:21781+5728？出于我的目的，我必须使用 EPSG:21781+5728 版本，那么我需要更改什么才能获得正确的像素值？