使用 GDAL 将立体投影转换为 epsg 28992

Question

我正在尝试使用 gdalwarp 将极地立体图像投影到墨卡托投影 epsg 28992（用于传单）

gdalinfo "HDF5:\"RAD_NL25_PCP_NA_202010271205.h5\"://image1/image_data"

Driver: HDF5Image/HDF5 Dataset
Files: RAD_NL25_PCP_NA_202010271205.h5
Size is 700, 765
Metadata:
  geographic_geo_column_offset=0
  geographic_geo_dim_pixel=KM,KM
  geographic_geo_number_columns=700
  geographic_geo_number_rows=765
  geographic_geo_par_pixel=X,Y
  geographic_geo_pixel_def=LU
  geographic_geo_pixel_size_x=1.0000035
  geographic_geo_pixel_size_y=-1.0000048
  geographic_geo_product_corners=0 49.362064 0 55.973602 10.856453 55.388973 9.0093002 48.895302
  geographic_geo_row_offset=3649.9819
  geographic_map_projection_projection_indication=Y
  geographic_map_projection_projection_name=STEREOGRAPHIC
  geographic_map_projection_projection_proj4_params=+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0
  image1_calibration_calibration_flag=Y
  image1_calibration_calibration_formulas=GEO = 0.500000 * PV + -32.000000
  image1_calibration_calibration_missing_data=0
  image1_calibration_calibration_out_of_image=255
  image1_image_bytes_per_pixel=1
  image1_image_geo_parameter=REFLECTIVITY_[DBZ]
  image1_image_product_name=RAD_NL25_PCP_H1.5_NA
  image1_image_size=535500
  image1_statistics_stat_max_value=44.5
  image1_statistics_stat_min_value=-31.5
  overview_hdftag_version_number=3.6
  overview_number_image_groups=1
  overview_number_radar_groups=3
  overview_number_satellite_groups=0
  overview_number_station_groups=0
  overview_products_missing=NA
  overview_product_datetime_end=27-OCT-2020;12:05:00.000
  overview_product_datetime_start=27-OCT-2020;12:05:00.000
  overview_product_group_name=RAD_NL25_PCP_NA
  radar1_radar_location=5.17834 52.101681
  radar1_radar_name=DeBilt
  radar1_radar_operational=0
  radar2_radar_location=4.7999701 52.953339
  radar2_radar_name=DenHelder
  radar2_radar_operational=1
  radar3_radar_location=5.1381001 51.836899
  radar3_radar_name=Herwijnen
  radar3_radar_operational=1
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  765.0)
Upper Right (  700.0,    0.0)
Lower Right (  700.0,  765.0)
Center      (  350.0,  382.5)
Band 1 Block=700x765 Type=Byte, ColorInterp=Undefined
  Metadata:
    image1_image_data_CLASS=IMAGE
    image1_image_data_PALETTE=
    image1_image_data_VERSION=1.2

gdalwarp -s_srs "+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0" -t_srs "epsg:28992" "HDF5:\"RAD_NL25_PCP_NA_202010271205.h5\"://image1/image_data" output.tif

导致此错误消息：

Processing HDF5:"RAD_NL25_PCP_NA_202010271205.h5"://image1/image_data [1/1] : 0ERROR 1: The transformation is already "north up" or a transformation between pixel/line and georeferenced coordinates cannot be computed for HDF5:"RAD_NL25_PCP_NA_202010271205.h5"://image1/image_data. There is no affine transformation and no GCPs. Specify transformation option SRC_METHOD=NO_GEOTRANSFORM to bypass this check.

我错过了什么？如果我设置标志 SRC_METHOD=NO_GEOTRANSFORM 我不会得到预期的结果

Answer 1

GDAL 中对 HDF5 数据的地理参考支持非常有限 (https://gdal.org/drivers/raster/hdf5.html)。因此，您的 HDF5 文件未被 GDAL 识别为地理参考数据。

我的解决方案是将您的命令分为两步：

使用 gdal_translate
使用 gdalwarp

首先，SRS的定义，地球的大小应该用m表示如下

+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378140 +b=6356750 +x_0=0 +y_0=0

接下来，假设元数据geographic_geo_product_corners表示图像四个角的经纬度，我们求出图像角的投影坐标

LL (0, -4415002.84084825)
UL (0, -3649999.11191775)
UR (700002.437056242, -3649999.05174429)
LR (700002.440031711, -4415003.15918867)

我用 gdaltransform 找到了这个。

gdaltransform -s_srs "+proj=latlong" -t_srs "+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378140 +b=6356750 +x_0=0 y_0=0"

所以，你问题中的转换可以通过以下命令完成

gdal_translate -of VRT -a_srs "+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378140 +b=6356750 +x_0=0 +y_0=0" -a_ullr 0 -3649999.05174429 700002.440031711 -4415003.15918867 "HDF5:\"RAD_NL25_PCP_NA_202010271205.h5\"://image1/image_data" /vsistdout/ | gdalwarp -of GTiff -t_srs epsg:28992 /vsistdin/ output.tif

有关如何将两个 gdal 命令与管道组合的信息，请参阅 How to convert projection of png tile from epsg:4326 to epsg:3857 by one command using gdal 中的答案。

使用 GDAL 将立体投影转换为 epsg 28992

Convert stereograhic projection to epsg 28992 with GDAL

gdal

proj

map-projections

mercator