1

我的 GeoServer 中有一个带有单通道（“灰色”）的高程栅格图层。“灰色”值是海拔值（有符号 int16）。我有 2 个客户：

  1. 第一个是按原样使用该高程值。
  2. 第二个期望得到[Mapbox Terrain-RGB 格式][1]

我不想将“灰度”格式转换为 Mapbox Terrain-RGB 格式并在 GeoServer 中保存重复数据。我正在考虑使用 SLD 样式和元素将高程值映射到适当的 RGB 值（在离散值之间进行梯度插值）。例如：

          <ColorMap>
            <ColorMapEntry color="#000000" quantity="-10000" />
            <ColorMapEntry color="#FFFFFF" quantity="1667721.5" />
          </ColorMap>

事实证明，上述示例并未涵盖所有颜色范围，而是仅创建灰度值。也就是说，它似乎对彼此独立的每种颜色（红色、绿色、蓝色）进行了插值。

知道如何让它插入这样的值：#000000、#000001、#000002、...、#0000FF、#000100、...、#0001FF、...、#FFFFFF？

德克萨斯州。[1]：https ://docs.mapbox.com/data/tilesets/reference/mapbox-terrain-rgb-v1/

4

2 回答 2

0

如果您拥有所需格式的高程数据，那么这是最简单的选择：它可以正常工作。但是，如果您想要一个更加定制化的解决方案，以下是我为使用 Mapbox Terrain-RGB 格式的项目所做的：

我有从深蓝色到浅蓝色到白色的颜色比例。我希望能够指定从浅蓝色到白色使用多少步（默认为 10）。此代码使用 GDAL Python 绑定。以下代码片段用于测试。它只是将颜色映射输出到 GeoTIFF 文件。要获取 0 到 1 之间的值，只需使用value *= 1/num_steps. 您可以在查找表中使用该值来获取 RGB 值。如果您只对输出颜色感兴趣，则可以忽略所有涉及gdal_translate. 颜色将自动存储在单波段 GeoTIFF 中。如果您确实想重复使用这些颜色，请注意此版本会忽略 alpha 值（如果存在）。您可以使用gdal_translate添加这些。该代码片段也可以在我的 gist 中找到。

import numpy as np
import gdal
from osgeo import gdal, osr

def get_color_map(num_steps):
    colors = np.zeros((num_steps, 3), dtype=np.uint8)
    colors[:, 0] = np.linspace(0, 255, num_steps, dtype=np.uint8)
    colors[:, 1] = colors[::-1, 0]
    return colors


ds = gdal.Open('/Users/myusername/Desktop/raster.tif')
band = ds.GetRasterBand(1) # Assuming single band raster
arr = band.ReadAsArray()
arr = arr.astype(np.float32)
arr *= 1/num_steps # Ensure values are between 0 and 1 (or use arr -= arr.min() / (arr.max() - arr.min()) to normalize to 0 to 1)


colors = get_color_map(num_steps) # Create color lookup table
colors[0] = [0, 0, 0] # Set black for no data so it doesn't show up as gray in the final product.


# Create new GeoTIFF with colors included (transparent alpha channel if possible). If you don't care about including the colors in the GeoTIFF, skip this.
cols = ds.RasterXSize
rows = ds.RasterYSize
out_ds = gdal.GetDriverByName('GTiff').Create('/Users/myusername/Desktop/raster_color.tif', cols, rows, 4)
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_ds.SetProjection(ds.GetProjection())
band = out_ds.GetRasterBand(1)
band.WriteArray(colors[arr]) # This can be removed if you don't care about including the colors in the GeoTIFF
band = out_ds.GetRasterBand(2)
band.WriteArray(colors[arr]) # This can be removed if you don't care about including the colors in the GeoTIFF
band = out_ds.GetRasterBand(3)
band.WriteArray(colors[arr]) # This can be removed if you don't care about including the colors in the GeoTIFF
band = out_ds.GetRasterBand(4)
alpha = np.zeros((rows, cols), dtype=np.uint8) # Create alpha channel to simulate transparency of no data pixels (assuming 0 is "no data" and non-zero is data). You can remove this if your elevation values are not 0.
alpha[arr == 0] = 255
band.WriteArray(alpha) # This can be removed if you don't care about including the colors in the GeoTIFF
out_ds.FlushCache()

使用具有多个值的调色板时，Rasterio 中也存在此问题。这是一个例子。

但是，如果您的栅格具有 n 维或者是掩码数组，则翻转操作可能会很棘手。这是基于此 stackoverflow 问题中的一个答案的解决方案：如何垂直翻转 2D NumPy 数组？.

于 2022-01-27T16:51:48.817 回答
0

我正在尝试做同样的事情，但没有运气，我认为它做不到......检查这个例子。这是一个带有 Mapbox 颜色编码的“渐变” [-10000, -5000, -1000, -500 ... 100000000000000000, 5000000000000000000, 1000000000000000000]。颜色进展/插值不是线性的，所以我认为它不能在 SLD 中模拟。

于 2022-01-29T06:01:52.030 回答