使用光栅重新投影坐标从规则网格创建不规则网格

Question

我已将 geotiff 文件加载到 xarray with a crs = EPSG:31467. I want to transform/reproject (don't know if there is a difference) these files into EPSG:4326. To do that, I use rasterio.warp.transform function which needs 1D arrays for x,y. To generate these i use numpy.meshgrid 并展开函数。这是我的数据的一个小例子：

import numpy

#Longitude and Latitude in EPSG:31467
lon = [3280914, 3281914, 3282914]
lat = [6103001, 6102001, 6101001]

#create 2d meshgrid
xv, yv = np.meshgrid(lon, lat)


xv, yv

    (array([[3280914, 3281914, 3282914],
        [3280914, 3281914, 3282914],
        [3280914, 3281914, 3282914]]),
 array([[6103001, 6103001, 6103001],
        [6102001, 6102001, 6102001],

        [6101001, 6101001, 6101001]]))

现在我有一个不同经度的序列 [3280914, 3281914, 3282914] 对于相同的纬度 [6103001, 6103001, 6103001] 当我现在使用 rasterio.transform(src_crs, dst_crs, x, y) 时，这些序列消失了，我不明白为什么？！

from rasterio.warp import transform

# Compute the lon/lat coordinates with rasterio.warp.transform
lon, lat = transform('EPSG:31467','EPSG:4326',
                     xv.flatten(), yv.flatten())



np.asarray(lon).reshape(3,3), np.asarray(lat).reshape(3,3)

> (array([[5.57397386, 5.58957607, 5.6051787 ],
>         [5.57473921, 5.59033795, 5.60593711],
>         [5.57550412, 5.5910994 , 5.60669509]]),  array([[55.00756605, 55.00800488, 55.00844171],
>         [54.9985994 , 54.99903809, 54.99947477],
>         [54.98963274, 54.99007128, 54.99050782]]))

np.unique(xv).shape, np.unique(yv).shape

> ((3,), (3,))

np.unique(lon).shape, np.unique(lat).shape

> ((9,), (9,))

要将重新映射的坐标改回 xarray，我必须在平等意义上获得相同的形状。哪个过程没看懂，是transform的功能还是projections的概念？

Answer 1

我不明白你在 np.asarray(lon).reshape(3,3)

之后到底想做什么

Which process I don't understand, is it the function of transform or the concept of projections?

看来你两个都不明白。

EPSG:31467 和 EPSG:4326 是根本上 不同类型的数据。 EPSG:31467实际上是纬向投影中的平面直角坐标系。 EPSG:4326根本就不是投影，是WGS-84大地坐标系中WGS-84椭球体的纯大地坐标。这里非常重要的是 EPSG:31467 中的相同坐标不必在 EPSG:4326 中相同。因为在 4326 中你的坐标是一个角度，而在 31467 中你的坐标是到赤道或假子午线的距离。这些系统中的轴不共线，并且与 子午线收敛 参数相关。因此，如果您在 31467 中更改 Norting 或 Easting，则纬度和对数都会发生变化。在这里你可以注意到蓝线（一个单元格是 31467 模拟）和黑线（整个网格是 4326 模拟）之间的角度 https://ru.wikipedia.org/wiki/%D0%A4%D0%B0%D0%B9%D0%BB:Soviet_topographic_map_kilometer_grid.svg

很容易检查转换是否正确 - 只需向后进行即可。

lon, lat = transform('EPSG:31467','EPSG:4326',
                     xv.flatten(), yv.flatten())
x_check, y_check = transform('EPSG:4326', 'EPSG:31467', lon, lat)

#we'll have some troubles because of computational errors, so let's round
x_check = [int(round(i, 0)) for i in x_check]

print(lon)
print(x_check)
print(xv.flatten())

>[5.574033001416839, 5.5896346633743175, 5.605236748547687, 5.574797816145165, 5.5903960110246524, 5.605994628800234, 5.5755622060626155, 5.591156935778857, 5.6067520880717225]
>[3280914, 3281914, 3282914, 3280914, 3281914, 3282914, 3280914, 3281914, 3282914]
>[3280914 3281914 3282914 3280914 3281914 3282914 3280914 3281914 3282914]

输出示例 transform() return 正是您所期望的 return。

下一个代码也能按预期工作（您可以将输出与上面的代码匹配）：

print(np.asarray(lon).reshape(3,3))
print(xv)

>[[5.574033   5.58963466 5.60523675]
> [5.57479782 5.59039601 5.60599463]
> [5.57556221 5.59115694 5.60675209]]
>[[3280914 3281914 3282914]
> [3280914 3281914 3282914]
> [3280914 3281914 3282914]]

我从未使用过 rasterio，因此无法为您提供可行的解决方案。一些注意事项：

我不知道为什么你需要网格来进行栅格转换
Rasterio 文档很清楚，并为您提供了解决方案：https://rasterio.readthedocs.io/en/latest/topics/reproject.html#reprojecting-a-geotiff-dataset
您可以直接在 crs 之间转换光栅。如果不在栅格中，请尝试 osgeo.gdal (gdal.Warp(dst_file, src_file, srcSRS='EPSG:31467', dstSRS='EPSG:4326')
请注意 重投影 和 定义栅格投影 之间的区别。首先更改图像，其次更改元数据。为了正确地进行直接变换，您的 GeoTIFF 必须在元数据中具有有效的投影定义（匹配实际光栅投影）
如果您不是在开发独立的应用程序，而只是需要重新投影 2-3 个栅格，请使用 QGIS 并且无需编码即可完成。在编码之前尝试理解 QGIS 中 2-3 个示例的大地测量概念也很有帮助。把它当作游乐场
如果你不是在开发独立的应用程序，你可以在 QGIS python API 中解决你的自动化任务。您可以使用 UI 测试工作流程，然后从 python 脚本中批量调用一些 QGIS/GDAL 工具。更重要的是 - rasterio 和所有其他软件包对于安装在 QGIS 上都是有价值的 python。当然，除非您正在创建 QGIS 插件，否则这不是部署的好主意
在EPSG:31467中0.001的坐标值为1mm。所以更精确是没有用的。在 EPSG:4326 中，1 度约为 111.1 公里（或 111.3*cos(lat)）。所以，你可以计算出有用的精确值。之后超过 4-5 位的所有内容。也可能没用

使用光栅重新投影坐标从规则网格创建不规则网格

Reproject coordinates with rasterio creates a irregular grid from an regular grid

python

projection

geospatial

python-xarray

rasterio