Advanced Examples¶
Optimize Transformations¶
Here are a few tricks to try out if you want to optimize your transformations.
Repeated transformations¶
Added in version 2.1.0.
If you use the same transform, using the pyproj.transformer.Transformer
can help
optimize your transformations.
import numpy
from pyproj import Transformer, transform
transformer = Transformer.from_crs(2263, 4326)
x_coords = numpy.random.randint(80000, 120000)
y_coords = numpy.random.randint(200000, 250000)
Example with pyproj.transformer.transform()
:
transform(2263, 4326, x_coords, y_coords)
Results: 160 ms ± 3.68 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Example with pyproj.transformer.Transformer
:
transformer.transform(x_coords, y_coords)
Results: 6.32 µs ± 49.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
Transforming with the same projections¶
pyproj skips noop transformations.
Transformation Group¶
Added in version 2.3.0.
The pyproj.transformer.TransformerGroup
provides both available
transformations as well as missing transformations.
Helpful if you want to use an alternate transformation and have a good reason for it.
>>> from pyproj.transformer import TransformerGroup
>>> trans_group = TransformerGroup("EPSG:4326","EPSG:2964")
>>> trans_group
<TransformerGroup: best_available=True>
- transformers: 8
- unavailable_operations: 1
>>> trans_group.best_available
True
>>> trans_group.transformers[0].transform(66, -153)
(149661.2825058747, 5849322.174897663)
>>> trans_group.transformers[1].transform(66, -153)
(149672.928811047, 5849311.372139239)
>>> trans_group.transformers[2].transform(66, -153)
(149748.32734832275, 5849274.621409136)
Helpful if want to check that the best possible transformation exists. And if not, how to get the missing grid.
>>> from pyproj.transformer import TransformerGroup
>>> tg = TransformerGroup("EPSG:4326", "+proj=aea +lat_0=50 +lon_0=-154 +lat_1=55 +lat_2=65 +x_0=0 +y_0=0 +datum=NAD27 +no_defs +type=crs +units=m", always_xy=True)
UserWarning: Best transformation is not available due to missing Grid(short_name=ntv2_0.gsb, full_name=, package_name=proj-datumgrid-north-america, url=https://download.osgeo.org/proj/proj-datumgrid-north-america-latest.zip, direct_download=True, open_license=True, available=False)
f"{operation.grids[0]!r}"
>>> tg
<TransformerGroup: best_available=False>
- transformers: 37
- unavailable_operations: 41
>>> tg.transformers[0].description
'axis order change (2D) + Inverse of NAD27 to WGS 84 (3) + axis order change (2D) + unknown'
>>> tg.unavailable_operations[0].name
'Inverse of NAD27 to WGS 84 (33) + axis order change (2D) + unknown'
>>> tg.unavailable_operations[0].grids[0].url
'https://download.osgeo.org/proj/proj-datumgrid-north-america-latest.zip'
Area of Interest¶
Added in version 2.3.0.
Depending on the location of your transformation, using the area of interest may impact which transformation operation is selected in the transformation.
>>> from pyproj.transformer import Transformer, AreaOfInterest
>>> transformer = Transformer.from_crs("EPSG:4326", "EPSG:2694")
>>> transformer
<Concatenated Operation Transformer: pipeline>
Description: Inverse of Pulkovo 1995 to WGS 84 (2) + 3-degree Gauss-Kruger zone 60
Area of Use:
- name: Russia
- bounds: (18.92, 39.87, -168.97, 85.2)
>>> transformer = Transformer.from_crs(
... "EPSG:4326",
... "EPSG:2694",
... area_of_interest=AreaOfInterest(-136.46, 49.0, -60.72, 83.17),
... )
>>> transformer
<Concatenated Operation Transformer: pipeline>
Description: Inverse of NAD27 to WGS 84 (13) + Alaska Albers
Area of Use:
- name: Canada - NWT; Nunavut; Saskatchewan
- bounds: (-136.46, 49.0, -60.72, 83.17)
Promote CRS to 3D¶
Added in version 3.1.
In PROJ 6+ you need to explicitly change your CRS to 3D if you have 2D CRS and you want the ellipsoidal height taken into account.
>>> from pyproj import CRS, Transformer
>>> transformer = Transformer.from_crs("EPSG:4326", "EPSG:2056", always_xy=True)
>>> transformer.transform(8.37909, 47.01987, 1000)
(2671499.8913080636, 1208075.1135782297, 1000.0)
>>> transformer_3d = Transformer.from_crs(
... CRS("EPSG:4326").to_3d(),
... CRS("EPSG:2056").to_3d(),
... always_xy=True,
...)
>>> transformer_3d.transform(8.37909, 47.01987, 1000)
(2671499.8913080636, 1208075.1135782297, 951.4265527743846)
Demote CRS to 2D¶
Added in version 3.6.
With the need for explicit 3D CRS since PROJ 6+, one might need to retrieve their 2D version, for example to create another 3D CRS compound between a 2D CRS and a vertical CRS.
>>> from pyproj import CRS, Transformer
>>> from pyproj.crs import CompoundCRS
>>> src_crs = CRS("EPSG:4979") # Any 3D CRS, here the 3D WGS 84
>>> vert_crs = CRS("EPSG:5773") # Any vertical CRS, here the EGM96 geoid
>>> dst_crs = CompoundCRS(src_crs.name + vert_crs.name, components=[src_crs.to_2d(), vert_crs])
>>> transformer_3d = Transformer.from_crs(src_crs, dst_crs, always_xy=True)
>>> transformer_3d.transform(8.37909, 47.01987, 1000)
(8.37909, 47.01987, 951.7851086745321)
Projected CRS Bounds¶
Added in version 3.1.
The boundary of the CRS is given in geographic coordinates. This is the recommended method for calculating the projected bounds.
>>> from pyproj import CRS, Transformer
>>> crs = CRS("EPSG:3857")
>>> transformer = Transformer.from_crs(crs.geodetic_crs, crs, always_xy=True)
>>> transformer.transform_bounds(*crs.area_of_use.bounds)
(-20037508.342789244, -20048966.104014594, 20037508.342789244, 20048966.104014594)
Multithreading¶
As of version 3.1, these objects are thread-safe:
If you have pyproj<3.1, you will need to create the object within the thread that uses it.
Here is a simple demonstration:
import concurrent.futures
from pyproj import Transformer
def transform_point(point):
transformer = Transformer.from_crs(4326, 3857)
return transformer.transform(point, point * 2)
with concurrent.futures.ThreadPoolExecutor(max_workers=10) as executor:
for result in executor.map(transform_point, range(5)):
print(result)
Optimizing Single-Threaded Applications¶
If you have a single-threaded application that generates many objects, enabling the use of the global context can provide performance enhancements.
For information about using the global context, see: Global Context
Here is an example where enabling the global context can help:
import pyproj
codes = pyproj.get_codes("EPSG", pyproj.enums.PJType.PROJECTED_CRS, False)
crs_list = [pyproj.CRS.from_epsg(code) for code in codes]
Caching pyproj objects¶
If you are likely to re-create pyproj objects such as pyproj.transformer.Transformer
or pyproj.crs.CRS
, using a cache can help reduce the cost
of re-creating the objects.
Transformer¶
from functools import lru_cache
from pyproj import Transformer
TransformerFromCRS = lru_cache(Transformer.from_crs)
Transformer.from_crs(2263, 4326) # no cache
TransformerFromCRS(2263, 4326) # cache
Try it:
from timeit import timeit
timeit(
"CachedTransformer(2263, 4326)",
setup=(
"from pyproj import Transformer; "
"from functools import lru_cache; "
"CachedTransformer = lru_cache(Transformer.from_crs)"
),
number=1000000,
)
timeit(
"Transformer.from_crs(2263, 4326)",
setup=("from pyproj import Transformer"),
number=100,
)
Without the cache, it takes around 2 seconds to do 100 iterations. With the cache, it takes 0.1 seconds to do 1 million iterations.
CRS Example¶
from functools import lru_cache
from pyproj import CRS
CachedCRS = lru_cache(CRS)
crs = CRS(4326) # no cache
crs = CachedCRS(4326) # cache
Try it:
from timeit import timeit
timeit(
"CachedCRS(4326)",
setup=(
"from pyproj import CRS; "
"from functools import lru_cache; "
"CachedCRS = lru_cache(CRS)"
),
number=1000000,
)
timeit(
"CRS(4326)",
setup=("from pyproj import CRS"),
number=1000,
)
Without the cache, it takes around 1 seconds to do 1000 iterations. With the cache, it takes 0.1 seconds to do 1 million iterations.
Debugging Internal PROJ¶
Added in version 3.0.0.
To get more debugging information from the internal PROJ code:
Set the
PROJ_DEBUG
environment variable to the desired level.Activate logging in pyproj with the devel DEBUG:
More information available here: https://docs.python.org/3/howto/logging.html
Here are examples to get started.
Add handler to the pyproj logger:
import logging console_handler = logging.StreamHandler() formatter = logging.Formatter("%(levelname)s:%(message)s") console_handler.setFormatter(formatter) logger = logging.getLogger("pyproj") logger.addHandler(console_handler) logger.setLevel(logging.DEBUG)
Activate default logging config:
import logging logging.basicConfig(format="%(levelname)s:%(message)s", level=logging.DEBUG)