Category Archives: fiona

What does it look like if you move all countries onto the same location?

2 Replies

Note from 2019-12-07: This is neither efficient nor up-to-date with modern PROJ. Better don’t copy and paste but just take non-projection/-transformation parts if you need them…

Sunday pre-lunch Python fun: What does it look like if you move all countries onto the same location?

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

import fiona
from shapely.geometry import *
from fiona.crs import from_string
from fiona.transform import transform_geom
from shapely.affinity import translate
 
def get_biggest_polygon(multipolygon):
    assert isinstance(geometry, MultiPolygon)
 
    max_area = 0
    biggest_polygon = None
    for polygon in geometry:
        if polygon.area > max_area:
            max_area = polygon.area
            biggest_polygon = polygon
    return biggest_polygon
 
def project_locally(geometry, from_crs):
    """centered on the centroid of the geometry"""
    # ugly because i map/shape back and forth, maybe try shapely instead
    lat = geometry.centroid.y
    lon = geometry.centroid.x
 
    to_crs = from_string("+proj=aeqd  +R=6371000 +lat_0={lat} +lon_0={lon}".format(lat=lat, lon=lon))
    reprojected_geometry = transform_geom(
        from_crs, 
        to_crs, 
        mapping(geometry)
    )
 
    return shape(reprojected_geometry)
 
with fiona.open("ne_10m_admin_0_countries.shp") as countries:
    centered_polygons = []
 
    for country in countries:
        geometry = shape(country['geometry'])
 
        # only use the biggest part of each country, otherwise everything sucks
        if isinstance(geometry, MultiPolygon):
            polygon = get_biggest_polygon(geometry)
        else:
            polygon = geometry
 
        # project nicely
        polygon = project_locally(polygon, countries.crs)
 
        # centering on 0,0 is simply moving the geometry by MINUS its x/y
        dx = -polygon.centroid.x
        dy = -polygon.centroid.y
        translated_polygon = translate(polygon, dx, dy)
 
        centered_polygons.append(translated_polygon)
 
with open("/tmp/outfile.wkt", "w") as sink:
    for polygon in centered_polygons:
        sink.write(polygon.wkt+"\n")

Use this in any way you like but please share your creations and code as well. :)

Some rough explanation: For each country I check if it is a multipolygon and if so, use only its biggest “sub”-polygon in the next steps. I then project the WGS84 coordinates to an Azimuthal Equidistant projection centered on the centroid of the polygon. That new geometry gets shifted to sit on the origin of the system. I collect all those polygons and write them as plain WKT to a file. Styling was done in QGIS.

And for the smart folk, the same but without local projection: