# Converting projections on Shapely

Posted on **2023-12-24**, by Racum.

Sometimes when dealing with coordinates, you may need to work in meters instead of degrees, but, since we don’t live in a flat earth, the calculation is not trivial. Thankfully, there are tools to help with the conversions.

## Understanting projections

Maps are a lie! We live in a spheric-ish, bumpy and irregular planet, but we need to convey geographical information on flat surfaces, and those things are not totally compatible; a map in a piece of paper of screen always is going to have some level of distortion.

There are many map projections to pick, each one with their own advantages and disadvantages: some preserve area, some preserve shapes, some are officially adopted by some countries, etc.

### Oblique Mercator

For the use case of performing operations in meters and converting back into degrees, the best projection is the Oblique Mercator, this guarantees a very low distortion near its anchor point. You probably already heard about the Mercator projection, very popular on printed maps; so, Oblique Mercator is a variant where the anchor point is arbitrary, while the *"normal"* Mercador has its anchor pinned at *“Null Island”* (latitude 0, longitude 0).

## Python libraries

### Shapely

Shapely is a very versatile geometry library, it contains geometric objects (Point, LineString, Polygon, etc) and a huge list of methods and operations to handle them.

But notice that Shapely is a **geometry** library, not a **geography** one! That means it is just a mathematical library, operating on a cartesian plane without attaching any units. Please refer to its very comprehensive manual if you want to get deeper.

#### Instalation:

```
$ pip install shapely
```

### pyproj

To use Shapely in a geography level, you need to combine it with pyproj: this is an interface for the PROJ library, the industry standard for all map projection transformations.

#### Instalation:

Before install pyproj via `pip`

, make sure you have PROJ installed, here are the instructions.

```
$ pip install pyproj
```

## Conversion functions

The conversion is not necessarily “hard”, but because PROJ is very flexible in its configuration, it can be intimidating trying to understand its documentation without some background. To make things easier, please use this snippet below, the function `get_projections()`

returns a tuple of two functions: the first converts geodesic shapes into the shapes in the cartesian plane, and the the second do the opposite.

```
from functools import partial
from pyproj import Proj, Transformer
from shapely.geometry import Point
from shapely.ops import transform
def get_projections(anchor_point=None):
if not anchor_point: # If no anchor, set 0/0 (normal Mercator).
anchor_point = Point(0, 0)
crs = {
'proj': 'omerc',
'lat_0': anchor_point.y,
'lonc': anchor_point.x,
'alpha': 1e-6,
'k': 1,
'gamma': 0.0,
'units': 'm',
'ellps': 'WGS84',
'no_defs': True,
}
proj = Proj(crs, preserve_units=True)
proj_inverse = Transformer.from_proj(proj, 'EPSG:4326', always_xy=True).transform
return partial(transform, proj), partial(transform, proj_inverse)
```

And this is how you use it: the `anchor_point`

is the geodesic point that will become our 0/0 point in the cartesian plane:

```
>>> center_of_berlin = Point(13.409408, 52.520842)
>>> brandenburg_gate = Point(13.377701, 52.516262)
>>> geodesic_to_cartesian, cartesian_to_geodesic = get_projections(center_of_berlin)
>>> geodesic_to_cartesian(brandenburg_gate)
<POINT (-2152.435 -509.177)> # In meters, from the center of Berlin.
>>> cartesian_to_geodesic(geodesic_to_cartesian(brandenburg_gate))
<POINT (13.378 52.516)> # Back to degrees.
```

## Examples

### Distance between points

Steps:

- Pick the source point as the anchor for
`get_projections()`

; - Transform both points into the cartesian plane with
`geodesic_to_cartesian()`

; - Run the
`source.distance(destination)`

.

Check the documentation on Shapely's distance() method.

```
>>> center_of_berlin = Point(13.409408, 52.520842)
>>> brandenburg_gate = Point(13.377701, 52.516262)
>>> geodesic_to_cartesian, _ = get_projections(center_of_berlin)
>>> cartesian_center_of_berlin = geodesic_to_cartesian(center_of_berlin)
>>> cartesian_brandenburg_gate = geodesic_to_cartesian(brandenburg_gate)
>>> cartesian_center_of_berlin.distance(cartesian_brandenburg_gate)
2211.8404871748057 # In meters.
```

### Area of a polygon

Steps:

- Load a Shapely polygon;
- Pick the polygon’s centroid as the anchor for
`get_projections()`

; - Transform it into the cartesian plane with
`geodesic_to_cartesian()`

; - Run
`polygon.area`

.

Check the documentation on Shapely's area property.

You can use both the constructor for the Polygon object or load it from GeoJSON with the shape() function:

```
from shapely.geometry import shape
museum_island = shape({
"type": "Polygon",
"coordinates": [
[
[13.393639, 52.522081],
[13.397824, 52.519285],
[13.403055, 52.512100],
[13.405072, 52.511486],
[13.411649, 52.513873],
[13.406418, 52.514692],
[13.401411, 52.519580],
[13.398160, 52.521513],
[13.393639, 52.522081]
]
]
})
```

Here is the shape above rendered in a map:

Please notice that if you get the `area`

from a geodesic object, the value will be in degrees, that makes no sense, since the atual area varies a lot if the shape is near the equator or near the poles.

```
>>> geodesic_to_cartesian, _ = get_projections(museum_island.centroid)
>>> museum_island.area
5.0027242999976146e-05 # ❌ In degrees, this is NOT useful!
>>> geodesic_to_cartesian(museum_island).area
377902.9629890913 # In square meters.
>>> geodesic_to_cartesian(museum_island).area / 1_000_000
0.3779029629890913 # In square kilometers.
```

### Radius polygon around a point

Steps:

- Pick the source point as the anchor or
`get_projections()`

; - Transform it into the cartesian plane with
`geodesic_to_cartesian()`

; - Apply a
`buffer()`

over it with your radius in meters (this returns a circular polygon); - Convert it back with
`cartesian_to_geodesic()`

; - (Optional): Use
`mapping()`

to convert it into GeoJSON.

Check the documentation on Shapely's buffer() method.

```
center_of_berlin = Point(13.409408, 52.520842)
geodesic_to_cartesian, cartesian_to_geodesic = get_projections(center_of_berlin)
cartesian_center_of_berlin = geodesic_to_cartesian(center_of_berlin)
cartesian_circle = cartesian_center_of_berlin.buffer(30_000, resolution=9) # 30Km, 36 vertices (4 * 9).
geodesic_circle = cartesian_to_geodesic(cartesian_circle)
```

This is how you convert the `geodesic_circle`

object into GeoJSON with the mapping() function:

```
import json
from shapely.geometry import mapping
circle_geojson = json.dumps(mapping(geodesic_circle))
```

The content of `circle_geojson`

looks like:

```
{
"type": "Polygon",
"coordinates": [
[
[13.851369517452731, 52.52001684395894], [13.84419333237743, 52.47322786545606],
[13.82384962144716, 52.427907468142855], [13.790992602415423, 52.385426101649315],
[13.746647020478344, 52.34706464968061], [13.692173951866176, 52.31397671258963],
[13.629227922924512, 52.287154992445906], [13.559706799470275, 52.267402631446686],
[13.485695926818149, 52.255310200330925], [13.40940799233307, 52.251238877281935],
[13.333120058076762, 52.255310201955346], [13.25910918610451, 52.267402634647084],
[13.189588063761999, 52.28715499712664], [13.126642036332331, 52.31397671861066],
[13.072168969589173, 52.34706465686153], [13.027823389824242, 52.38542610977443],
[12.994966373204546, 52.427907476967505], [12.974622664855207, 52.47322787471379],
[12.967446482452836, 52.52001685336939], [12.97369737930227, 52.5668552183338],
[12.993227388561863, 52.612317819578806], [13.025480427930086, 52.655017058880745],
[13.069504604188237, 52.69364577761025], [13.123977621128224, 52.727018141547795],
[13.187244976593817, 52.7541071425489], [13.25737006180486, 52.77407739795136],
[13.332194687865442, 52.786312051983195], [13.409408007761378, 52.790432768978405],
[13.48662132741714, 52.7863120503392], [13.561445952765052, 52.774077394714226],
[13.63157103681411, 52.754107137818636], [13.694838390706073, 52.72701813547045],
[13.749311405711442, 52.693645770373024], [13.79333557973581, 52.65501705070612],
[13.825588616641758, 52.61231781071741], [13.84511862328763, 52.5668552090565],
[13.851369517452731, 52.52001684395894]
]
]
}
```

And this is the same polygon rendered in a map:

## Conclusion

I hope you can find this snippet useful; feel free to use it on your projects.

The examples in this article are not even scratching the surface on Shapely, I recommend a back-to-back read on its documentation. You probably don’t need to dive as deep on PROJ, unless your project has a custom required standard to follow.

For things more advanced, the same patterns stands: **convert to cartesian** ➤ **perform operations** ➤ **convert to geodesic**.