Tes Kecepatan Proyeksi Geometri

Akhir-akhir ini saya telah menggunakan kelas proyeksi OGR yang datang dengan ogr / gdal, tetapi pyproj direkomendasikan kepada saya, jadi saya pikir saya akan mencobanya. Untuk membantu saya memutuskan apakah saya harus beralih, saya melakukan tes kecepatan. Berikut ini adalah contoh kecil (hampir) direproduksi yang saya buat untuk menguji keduanya. Saya tidak yakin apakah tes ini benar-benar adil, jadi komentar dan rekomendasi dipersilahkan!

Impor terlebih dahulu, untuk memastikan kami memulai dengan bidang permainan level:

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

Karena saya menyimpan geometri berbentuk dalam panda 'Seri', fungsi-fungsi perlu dikerjakan Series.apply(). Di sini saya mendefinisikan dua fungsi (satu menggunakan 'ogr.osr' dan satu menggunakan 'pyproj') untuk melakukan transformasi koordinat di dalam panggilan ke 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)

Masing-masing fungsi ini mengambil kode EPSG sebagai input untuk sistem referensi koordinat asal dan tujuan. Kedua perpustakaan menawarkan cara-cara alternatif untuk mendefinisikan informasi proyeksi, tetapi kode EPSG membuat kode tersebut cukup sederhana untuk saat ini.

Inilah hasilnya, menggunakan %timeitfungsi ajaib di 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

Jelas versi 'ogr.osr' lebih cepat, tetapi apakah ini perbandingan yang adil? Versi 'pyproj' dilakukan pada titik-titik individual, dan sebagian besar dijalankan dengan Python, sedangkan versi 'ogr.osr' beroperasi langsung pada objek geometri OGR. Apakah ada cara yang lebih baik untuk membandingkan ini?

Pyproj adalah ekstensi Python C yang menggunakan perpustakaan PROJ4 dan osgeo.ogr adalah ekstensi Python C yang menggunakan perpustakaan PROJ4. Jika Anda hanya mempertimbangkan proyeksi koordinat, dalam tes paling adil mereka akan hampir sama.

Transformasi Pyproj dapat beroperasi pada array nilai koordinat, jadi Anda hanya perlu memanggilnya sekali per baris atau dering alih-alih untuk setiap pasangan. Ini harus mempercepat beberapa hal. Contoh: https://gist.github.com/sgillies/3642564#file-2-py-L10 .

(Perbarui) Juga, Shapely menyediakan fungsi yang mengubah geometri di 1.2.16:

Help on function transform in module shapely.ops:

transform(func, geom)
    Applies `func` to all coordinates of `geom` and returns a new
    geometry of the same type from the transformed coordinates.

    `func` maps x, y, and optionally z to output xp, yp, zp. The input
    parameters may iterable types like lists or arrays or single values.
    The output shall be of the same type. Scalars in, scalars out.
    Lists in, lists out.

    For example, here is an identity function applicable to both types
    of input.

      def id_func(x, y, z=None):
          return tuple(filter(None, [x, y, z]))

      g2 = transform(id_func, g1)

    A partially applied transform function from pyproj satisfies the
    requirements for `func`.

      from functools import partial
      import pyproj

      project = partial(
          pyproj.transform,
          pyproj.Proj(init='espg:4326'),
          pyproj.Proj(init='epsg:26913'))

      g2 = transform(project, g1)

    Lambda expressions such as the one in

      g2 = transform(lambda x, y, z=None: (x+1.0, y+1.0), g1)

    also satisfy the requirements for `func`.