Más

¿Cómo calcular el cuadro delimitador en coordenadas proyectadas?

¿Cómo calcular el cuadro delimitador en coordenadas proyectadas?


Estoy usando GDAL en Python para reproyectar rásteres. Me gustaría poder tomar un ráster existente y proyectarlo en un nuevo ráster para mi proyección y tamaño de píxel de salida de elección.

Me preocupa cómo calcular el cuadro delimitador en coordenadas proyectadas en las que, a su vez, basaré la geotransformación.gdaltiene una función llamadagdal.AutoCreateWarpedVRTque aparentemente crea un ráster de memoria desde el que puedo leer directamente la geotransformación. Pero a veces tengo rásteres MUY grandes que ni siquiera caben en la memoria principal, odiaría crear un ráster basado en memoria solo para obtener la geotransformación.

¿Es suficiente transformar los vértices del cuadro delimitador original en el sistema de coordenadas proyectadas y luego de alguna manera convertir ese cuadro en un cuadrado? (Estoy seguro de que la caja estará distorsionada).


Una forma segura de obtener la nueva extensión exacta es

  1. crear un polígono a partir de la extensión del ráster con gdaltindex
  2. densificar la geometría del polígono
  3. reproyecte el polígono al CRS de destino con ogr2ogr
  4. calcular la extensión en el nuevo CRS

Han pasado unos meses desde que publiqué esto. @ michael-miles-stimson sugirió que la forma correcta de hacer esto sería muestrear puntos a lo largo de cada borde del cuadro delimitador, proyectarlos y luego construir un cuadro delimitador abarcador en esos puntos. El código desplegable de Python / GDAL a continuación hace esto y le permite seleccionar la cantidad de puntos de interpolación que le gustaría muestrear a lo largo del borde.

import numpy from osgeo import osr def transform_bounding_box (bounding_box, base_epsg, new_epsg, edge_samples = 11): "" "Transforma el cuadro delimitador de entrada en proyección de salida. Esta transformación explica el hecho de que el cuadro delimitador cuadrado reproyectado podría deformarse en la nueva coordenada sistema. Para tener en cuenta esto, la función muestrea puntos a lo largo de los bordes del cuadro delimitador original e intenta hacer el cuadro delimitador más grande alrededor de cualquier punto transformado en el borde, ya sean esquinas o bordes deformados. Parámetros: bounding_box (lista): una lista de 4 coordenadas en el sistema de coordenadas 'base_epsg' que describe el límite en el orden [xmin, ymin, xmax, ymax] base_epsg (int): el código EPSG del sistema de coordenadas de entrada new_epsg (int): el código EPSG del sistema de coordenadas de salida deseado edge_samples ( int): el número de puntos interpolados a lo largo de cada borde del cuadro delimitador para muestrear a lo largo. Un valor de 2 tomará una muestra solo de las esquinas, mientras que un valor de 3 también tomará una muestra de las esquinas y el punto medio. Devuelve: A lista de la forma [xmin, ymin, xmax, ymax] que describe el cuadro delimitador de ajuste más grande alrededor del cuadro delimitador deformado original en el sistema de coordenadas 'new_epsg'. "" "base_ref = osr.SpatialReference () base_ref. , cuadro_delimitador [3])) p_1 = numpy.array ((cuadro_delimitador [0], cuadro_delimitador [1])) p_2 = numpy.array ((cuadro_delimitador [2], cuadro_delimitador [1])) p_3 = numpy.array (( bounding_box [2], bounding_box [3])) def _transform_point (punto): trans_x, trans_y, _ = (transformer.TransformPoint (* point)) return (trans_x, trans_y) # Esta lista de comprensión itera sobre cada borde del cuadro delimitador , # divide cada borde en 'edge_samples' número de puntos, luego reduce # esa lista a un 'límite_fn' apropiado dado el borde. # Por ejemplo, el borde izquierdo debe ser la coordenada x mínima, por lo que # generamos el número de 'edge_samples' de puntos entre el punto superior izquierdo y el # punto inferior izquierdo, transfórmalos todos al nuevo sistema de coordenadas # y luego obtén la coordenada x mínima "min (p [0]…)" de la b atch. transform_bounding_box = [bounding_fn ([_transform_point (p_a * v + p_b * (1 - v)) para v en numpy.linspace (0, 1, edge_samples)]) para p_a, p_b, bounding_fn en [(p_0, p_1, lambda p_list : min ([p [0] para p en p_list])), (p_1, p_2, lambda p_list: min ([p [1] para p en p_list])), (p_2, p_3, lambda p_list: max ([ p [0] para p en p_list])), (p_3, p_0, lambda p_list: max ([p [1] para p en p_list]))]] return transform_bounding_box