¿Cómo añadir el sistema de coordenadas a los archivos de máscara de imágenes Spot?

Question

Trabajo con imágenes Spot 6 y 7 suministradas en formato DIMAP. He descubierto cómo extraer las bandas multiespectrales y pancromáticas en geotiff con gdal-translate para que sea más fácil trabajar con ellos.

En cada archivo hay una serie de archivos de máscara en gml que puedo leer con ogrinfo y Qgis. Sin embargo, los archivos de máscara no tienen un sistema de coordenadas, por lo que no puedo utilizarlos con las imágenes.

Según el informe ogrinfo, parece que el GML utiliza las dimensiones en píxeles de las filas y columnas de la imagen. (La imagen de origen que coincide es de 9652 x 57083).

$ ogrinfo \SPOT6_sample_roi.gml maskfeature
INFO: Open of `\SPOT6_sample_roi.gml'
      using driver `GML' successful.
Metadata:
  NAME=Area of interest mask for product id SPOT6_MS_201308032015087_SEN_SPOT6_20160316_1601281mdxzlrvssw12_1

Layer name: MaskFeature
Geometry: Polygon
Feature Count: 1
Extent: (1.000000, 1.000000) - (9653.000000, 57084.000000)
Layer SRS WKT:
(unknown)
gml_id: String (0.0) NOT NULL
maskType: String (18.0)
OGRFeature(MaskFeature):0
  gml_id (String) = REGION_OF_INTEREST-0
  maskType (String) = REGION_OF_INTEREST
  POLYGON ((9645.1767578125 6.41328716278076,9645.162109375 5.32240867614746,9645.03125 4.30024194717407,9644.7841796875 3.2741334438324,9644.5390625 2.39344930648804,9644.181640625 1.72693908214569,9643.7666015625 1.20388793945312,9643.4931640625 1.0,8939.99609375 1.0,1.0 4.3671669960022,1.0 28542.5,1 57084,9653 57084,9653.0 45960.265625,9645.1767578125 6.41328716278076))

El informe gdalinfo de la imagen de origen también muestra las coordenadas de los píxeles para el sistema de coordenadas:

Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,57083.0)
Upper Right ( 9652.0,    0.0)
Lower Right ( 9652.0,57083.0)
Center      ( 4826.0,28541.5)

Sin embargo, también tiene metadatos RPC que parecen ser suficientes para que Qgis y ArcMap/Pro lo muestren en la ubicación geográfica correcta:

RPC Metadata:
  HEIGHT_OFF=500.0
  HEIGHT_SCALE=500.0
  LAT_OFF=64.90742355
  LAT_SCALE=1.71583845
  ...snip...
  SAMP_OFF=4825
  SAMP_SCALE=4826.0

Archivos de muestra en https://drive.google.com/drive/folders/119QEECJ42FKt0A9mq55rRmGsT9-5nSZq?usp=sharing (La imagen ha sido redimensionada al 10% de su tamaño original).

¿Cómo puedo combinar la información del sistema de coordenadas rasterizado con los archivos de máscara para poder utilizarlos juntos?

Answer 1

Un pequeño avance: podemos utilizar gdaltransform -rpc para obtener resultados lat-long (o cualquier otro sistema de coordenadas) para las ubicaciones de los píxeles de la imagen:

$ gdaltransform -rpc DS_SPOT6_201308032015087_LM1_FR1_FR1_FR1_W137N65_10315\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A\DIM_SPOT6_MS_201308032015087_SEN_1.XML
Enter column line values separated by space, and press Return.
0.0 57083.0
-137.766376889524 63.1915449659547 0
9652.0    0.0
-136.420879379204 66.6213864498263 0

Estas son las coordenadas de las esquinas del informe gdalinfo, y coinciden con la imagen tal y como la muestra ArcMap.

Así que en python podríamos hacer un bucle a través de los vértices de cada par de coordenadas del polígono GML y transformar usando el método anterior, pero parece que debería haber un camino más directo.