Últimamente he estado usando las clases de proyección OGR que vienen con ogr / gdal, pero me recomendaron pyproj, así que pensé en probarlo. Para ayudarme a decidir si debería hacer el cambio, hice una prueba de velocidad. El siguiente es un pequeño (casi) ejemplo reproducible que se me ocurrió para probar los dos. No estoy seguro de si esta prueba es totalmente justa, ¡así que los comentarios y recomendaciones son bienvenidos!
Importa primero, para asegurarnos de comenzar con un campo de juego nivelado:
from pandas import Series # This is what our geometries are stored in
from shapely import wkb
import functools
from shapely.geometry import shape
from osgeo import ogr
# The following two imports are the important ones
from pyproj import Proj, transform
from osgeo.osr import SpatialReference, CoordinateTransformation
Debido a que estoy almacenando las geometrías bien formadas en una 'Serie' de pandas, las funciones deben funcionar Series.apply()
. Aquí defino dos funciones (una usando 'ogr.osr' y otra usando 'pyproj') para realizar transformaciones de coordenadas dentro de una llamada a Series.apply()
:
def transform_osr(geoms, origin, target):
target_crs = SpatialReference()
target_crs.ImportFromEPSG(origin)
origin_crs = SpatialReference()
origin_crs.ImportFromEPSG(origin)
transformer = CoordinateTransformation(origin_crs, target_crs)
def trans(geom):
g = ogr.CreateGeometryFromWkb(geom.wkb)
if g.Transform(transformer):
raise Exception("Ahhhh!")
g = wkb.loads(g.ExportToWkb())
return g
return geoms.apply(trans)
def transform_pyproj(geoms, origin, target):
target = Proj(init="epsg:%s" % target)
origin = Proj(init="epsg:%s" % origin)
transformer = functools.partial(transform, origin, target)
def trans(geom):
interface = geom.__geo_interface__
geom_type = interface['type']
coords = interface['coordinates']
result = apply_to_coord_pairs(transformer, coords)
return shape({'coordinates':result, 'type':geom_type})
def apply_to_coord_pairs(fun, geom):
return [not all([hasattr(y, "__iter__") for y in x]) and \
fun(*x) or apply_to_coord_pairs(fun, x) for x in geom]
return geoms.apply(trans)
Cada una de estas funciones toma un código EPSG como entrada para los sistemas de referencia de coordenadas de origen y destino. Ambas bibliotecas ofrecen formas alternativas de definir información de proyección, pero los códigos EPSG mantienen el código bastante simple por ahora.
Aquí están los resultados, usando la %timeit
función mágica en ipython:
In [1]: %timeit transform_pyproj(l, 29903, 4326)
1 loops, best of 3: 641 ms per loop
In [2]: %timeit transform_osr(l, 29903, 4326)
10 loops, best of 3: 27.4 ms per loop
Claramente, la versión 'ogr.osr' es más rápida, pero ¿es una comparación justa? La versión 'pyproj' se realiza en puntos individuales y se ejecuta principalmente en Python, mientras que la versión 'ogr.osr' opera directamente en el objeto de geometría OGR. ¿Hay una mejor manera de comparar esto?
fuente
projected_coords = numpy.vstack(pyproj.transform(origin, target, *numpy.array(geom).T)).T