0

我有一个 Geotiff，我在瓷砖地图上显示，但它稍微偏南。例如，在这个屏幕截图中，图像的边缘应该是国家边界所在的位置，但它有点偏南：

问题截图

这是代码的相关部分：

tiff_rio_500 = rioxarray.open_rasterio('/content/mw/mw_dist_to_light_at_all_from_light_mask_mw_cut_s3_500.tif')
dataarray_500 = tiff_rio_500[0]
dataarray_500_meters = dataarray_500.copy()
dataarray_500_meters['x'], dataarray_500_meters['y'] = ds.utils.lnglat_to_meters(dataarray_500.x, dataarray_500.y)
hv_dataset_500_meters = hv.Dataset(dataarray_500_meters, name='nightlights', vdims='cumulative_cost')
hv_tiles_osm_bokeh = hv.element.tiles.OSM().opts(width=1000, height=800)
hv_image_500_meters_bokeh = hv.Image(hv_dataset_500_meters, kdims=['x', 'y'], vdims=['cumulative_cost'], rtol=1).opts(cmap='inferno_r')
hv_combined_osm_500_meters_bokeh = hv_tiles_osm_bokeh * hv_image_500_meters_bokeh
hv_combined_osm_500_meters_bokeh

你可以在 google colab 上看到 live notebook。

现在，当人们不将地图转换为 Web 墨卡托时，这不是常见的“一切都结束了”的问题。它几乎是完美的，但事实并非如此。

Geotiff 是一个地球引擎出口。这是它在 Earth Engine 中最初的样子（参见实时代码）：
Earth Engine 中一切正常的屏幕截图
正如您所见，图像随处可见。

起初，我怀疑可能是导出出错了，或者谷歌地图瓦片集有些不同，但不是，如果我在 Windows 笔记本电脑上的 QGis 应用程序中打开相同的导出 Tiff 并在同一个 OSM 瓦片地图上查看它在 colab 笔记本中，它看起来不错：
QGis中一切正常的截图
好的，图像没有完美地遵循边界，但我知道为什么而且这不相关（我过度简化了国家边界几何）。关键是，它被投影到正确的位置。因此，基于此，tiff 包含正确的信息，它可以显示在与 OSM tilemap 中的边界相同的位置，但在我的 Holoviews-Datashader-Bokeh 项目中，它仍然略有偏差。

知道为什么会这样吗？

4

1 回答 1

0

我从一位开发人员那里得到了关于 Holoviz Discourse 的答案。看到推荐的功能实际上是如何无证的，我在这里复制它，以防有人寻找一种简单的方法来加载 geotiff 并添加到 Holoviews/Geoviews 中的 tilemap：

https://discourse.holoviz.org/t/geotiff-overlay-position-is-slightly-off-on-holoviews-bokeh-tilemap/2071

philippjfr
我不希望手动转换坐标工作得特别好。虽然对于精确的坐标变换来说，它对权重的依赖性要大得多，但我建议使用GeoViews。

img = gv.util.load_tiff('/content/mw/mw_dist_to_light_at_all_from_light_mask_mw_cut_s3_500.tif') gv.tile_sources.OSM() * img.opts(cmap='inferno_r')

编辑：现在可能不想使用 Geoviews，因为它有一个非常重的依赖链，需要很大的耐心和运气才能正确设置它。幸运的是 rioxarray（通过 rasterio）有一个重新投影的工具，只需将“.rio.reproject('EPSG:3857')”附加到第一行，然后您就不必使用不用于此目的的 lnglat_to_meters。

所以修正后的代码变成：

tiff_rio_500 = rioxarray.open_rasterio('/content/mw/mw_dist_to_light_at_all_from_light_mask_mw_cut_s3_500.tif').rio.reproject('EPSG:3857')
hv_dataset_500_meters = hv.Dataset(tiff_rio_500[0], name='nightlights', vdims='cumulative_cost')
hv_tiles_osm_bokeh = hv.element.tiles.OSM().opts(width=1000, height=800)
hv_image_500_meters_bokeh = hv.Image(hv_dataset_500_meters, kdims=['x', 'y'], vdims=['cumulative_cost'], rtol=1).opts(cmap='inferno_r')
hv_combined_osm_500_meters_bokeh = hv_tiles_osm_bokeh * hv_image_500_meters_bokeh
hv_combined_osm_500_meters_bokeh

现在与 Geoviews 解决方案（据说会自动处理所有内容）相比，此解决方案有一个缺点，即如果您使用悬停工具提示在鼠标光标下显示值和坐标，坐标将显示在新投影的 web 墨卡托系统中数百万米，而不是预期的度数。解决方案超出了此答案的范围，但我刚刚完成了一个详细的分步指南，其中也包含了一个解决方案，一旦发布，我将在此处链接。当然，如果您不使用 Hover Tooltip，上面的代码将非常适合您，无需任何修改。

于 2021-04-08T14:14:10.430 回答