**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:

Pingback: GeoArt! | geoObserver

Pingback: 2 x GeoArt! | geoObserver