¿Cómo convierto coordenadas afines a lat / lng?

14

Soy extremadamente nuevo en SIG.

Estoy usando gdalpara leer en un mapa de uso / cobertura del suelo y necesito seleccionar el lat / lng de ciertos tipos de cobertura del suelo para indexar en un conjunto de datos diferente que se expresa solo en lat / lng. Desafortunadamente, no entiendo la forma de las coordenadas X e Y que me ha dado de la transformación geográfica, en concreto el originXy originYa continuación:

geotransform = dataset.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]

Imprimir estos valores me da coordenadas como (447466.693808, 4952570.40529). ¿Cómo se relacionan con la latitud y longitud originales?

Editar:

Aquí hay un ejemplo simple de Python que me dio lo que estaba buscando:

srs = osr.SpatialReference()
srs.ImportFromWkt(dataset.GetProjection())

srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srs,srsLatLong)
print ct.TransformPoint(originX,originY)

Robado de: tolatlong.py

Rico
fuente
Parece que sus datos están proyectados (por ejemplo, UTM) y necesitará saber cuál es la proyección para "desproyectarla" de nuevo a coordenadas largas / largas.
@Dan Gracias, así que sé que puedo obtener la proyección a través de ay dataset.GetProjectionRef()descubrir que estoy usando "UTM Zone 10", pero ¿y luego qué? Estoy buscando métodos como "desproyectar", pero estoy llegando a nulo.
Rico
perdón por el uso del término unproject (entre comillas) ya que los datos en grados decimales no se proyectan, pero si desea recuperar los datos proyectados a grados decimales de cualquier proyección dada, entonces debe (tenga en cuenta las comillas) "proyecto" vuelve a un sistema de coordenadas geográficas, también conocido como datos de grados decimales.
2
Este hilo (más reciente) proporciona un ejemplo explícito y otra solución: gis.stackexchange.com/questions/8430/…
whuber

Respuestas:

10

gdal_translate reproyectará sus datos de cualquier proyección en cualquier otra (en este caso desea EPSG: 4326) usando:

gdal_translate -a_srs epsg:4326 srcfile new_file 

o podría usar gdaltrasform para convertir los puntos (y estoy seguro de que también puede acceder a eso desde Python (?))

Ian Turton
fuente
(+1) No sé por qué esto fue rechazado, Ian, porque me parece una especie de operación que será necesaria después de que se aplique la transformación afín.
whuber
gdal_translate no funcionará como lo tiene (o más bien simplemente asignaría EPSG: 4326 a los datos), debe deformar los datos con algo como: gdalwarp -s_srs EPSG: 32610 -t_srs EPSG: 4326 src dest Si el los datos ya tienen metadatos de proyecto, puede deshacerse del parámetro -s_srs.
MerseyViking
Tal vez esto es lo que busco. ¿Es el "epsg: 4326" una transformación a lat / lng global? Creo que veo algunas clases donde puedo proporcionar una proyección, estaba pensando "WGS: 84", pero no sé la diferencia.
Rico
1
@ Rich Haga clic en la etiqueta "sistema de coordenadas" en su pregunta, ordene la página resultante por votos y comience a leer. Muchas de sus preguntas serán respondidas dentro de unos minutos.
whuber
7

La geotransformación se documenta en https://gdal.org/user/raster_data_model.html . La idea es que tome coordenadas (x, y) directamente del conjunto de datos, aplique una transformación lineal para obtener (u, v) con

u = a*x + b*y
v = c*x + d*y

(puede tomar esto como la definición de una transformación lineal), luego cambie el origen agregando geotransform [0] a u y geotransform [3] a v. Eso da la "transformación afín" de (x, y). Realmente tiene la intención de rotar, cambiar la escala, quizás corregir un poco para algunos errores de sesgo y reubicar las coordenadas específicas de los datos (x, y) para que coincidan con un sistema de coordenadas conocido. Se supone que el resultado produce coordenadas proyectadas. Esto simplemente significa que hay un procedimiento matemático que toma (longitud, latitud) y los convierte en las coordenadas conocidas: esto se llama "proyección". "Desproyectar" es hacer lo contrario; entonces, si sabes qué proyección es necesaria, la aplicas a las coordenadas transformadas afines (x, y) para obtener la latitud y longitud.

Por cierto, los valores de las constantes a, b, c, d están dados por las entradas 1, 2, 4 y 5 en la matriz de geotransformación.

whuber
fuente
1
Leí la sección de geotransformación del enlace que proporcionó, pero no creo que esto sea lo que busco. Lo siento si estoy siendo obtuso, pero ¿no describe esto cómo proyectar los índices x e y (0 a XSize o YSize) a las coordenadas proyectadas? Tengo curiosidad por cómo traduce las coordenadas proyectadas en latitud y longitud. Como ejemplo simple, la geotransformada define el origen x / y de la esquina superior izquierda del ráster en las coordenadas proyectadas (447466.693808, 4952570.40529). ¿Cómo volvería a convertir esas coordenadas en lat / lng (algo así como lat: 44, lng: -122)?
Rico
@ Rich No, ese enlace no describe una proyección. Solo describe un cambio de coordenadas entre dos mapas, no entre el mundo y un mapa. La no proyección convierte la transformación afín de las coordenadas en (lat, lon). Tú lo haces no quiere escribir el código para que si se puede evitar! La no proyección es casi siempre una cuestión de identificar qué proyección se necesita y llamar al software adecuado para que haga el trabajo por usted (como se sugiere en la respuesta de @iant).
whuber
El enlace está roto. @whuber, ¿tengo razón en que la transformación afín (de GetGeoTransform de GDAL) es una transformación entre píxeles y coordenadas CRS? Entonces, por ejemplo, si los datos están en una proyección GDA, los puntos deben transformarse de WGS84 / lat / lon a GDAXX usando, por ejemplo, CoordianteTransform, y luego a píxeles usando la afín GeoTransform. Este proceso de dos pasos es algo que me ha estado tropezando bastante, y sospecho que también está causando problemas con el OP. No ayudó que xarray / rasterio parece tener un error que arruina los valores de GeoTransform: GeoTransform github.com/pydata/xarray/issues/3185
naught101
@naught Encontré la nueva URL y actualicé mi publicación.
whuber
0

Puedes usar lo siguiente:

import gdal, numpy as n
def coord(file):
    padfTransform = file.GetGeoTransform()
    indices = n.indices(file.ReadAsArray().shape)
    xp = padfTransform[0] + indices[1]*padfTransform[1] + indices[1]*padfTransform[2]   
    yp = padfTransform[3] + indices[0]*padfTransform[4] + indices[0]*padfTransform[5]  
    return xp,yp
file = gdal.Open('Yourfile.tif') # A GeoTiff file
x,y = coord(file)

coord devolverá la longitud (x) y la latitud (y) de todos los píxeles. Tenga en cuenta que las coordenadas son de la esquina izquierda de un píxel

Ishan Tomar
fuente
¿No deberían ser índices [1] * padfTransform [1] + índices [0] * padfTransform [2] y viceversa para el yp?
CMCDragonkai