Notebooks from applied-cs/data-science@feb99093

4d2059da · Christof Kaufmann · 71b153b1 · 4d2059da
Commit 4d2059da authored 11 months ago by Christof Kaufmann
--- a/10-geodaten/01-projections.ipynb
+++ b/10-geodaten/01-projections.ipynb
@@ -86,7 +86,7 @@
    "\n",
    "if scaledown:\n",
    "    # split MultiPolygons into separate Polygons, e. g. 1 MultiPolygon of Germany → 6 Polygons (6 rows with the same data)\n",
-    "    countries = countries.explode()\n",
+    "    countries = countries.explode(index_parts=False)\n",
    "\n",
    "    # scale down each poly\n",
    "    scale = 1 / np.sqrt(ratio)\n",
@@ -112,7 +112,7 @@
    "    attribution=countries.crs.name\n",
    ")"
   ],
-   "id": "0005-87f4bcf5159c2d37428cabbcd55f47bf4cd95bbb86ed15eced56f1b465b"
+   "id": "0005-5b4adbf9e675eff3ed4415e452a13cfb859e618ba5637ea6a5290655f08"
  },
  {
   "cell_type": "markdown",

 %% Cell type:markdown id:0002-a8f53080980ed874d387cc3bfb8c5d5ed65d7fdb4139eb4ed15105d383c tags:

 # Flächenverzerrung bei Projektionen

 In dieser Demo kann man sich die Flächenverzerrungen bei verschiedenen
 Projektionen auf einer interaktiven Karte anschauen.

 Zunächst laden wir eine GeoJSON-Datei mit groben Ländergrenzen herunter
 und speichern sie lokal. Und wir laden auch eine mehr oder weniger gute
 Datei mit den Referenzgrößen. Mit etwas Vorverarbeitung können wir die
 Datensätze in `countries` zusammenführen, was wir aber erst im nächsten
 Schritt machen.

 %% Cell type:code id:0003-1224f13613ae889be0984402ff1578d533b9774aae88ca6305a774f6860 tags:

 ``` 
 import numpy as np
 import os
 import geopandas as gpd
 from shapely import affinity
 from sklearn.datasets import fetch_openml

 path = 'geo_countries.geojson'
 if not os.path.exists(path):
    url = 'https://d2ad6b4ur7yvpq.cloudfront.net/naturalearth-3.3.0/ne_50m_admin_0_countries.geojson'
    geo_countries = gpd.read_file(url)
    geo_countries.to_file(path)      # save local copy

 geo_countries = gpd.read_file(path)  # load from local copy
 geo_countries = geo_countries[['name', 'iso_a3', 'geometry']].set_index('iso_a3')

 countries_df, _ = fetch_openml('Countries-of-the-World', as_frame=True, return_X_y=True)
 countries_df = countries_df.rename(columns={'Area_(sq._mi.)': 'ref_area', 'Country': 'name'})  # unit was km² already
 countries_df = countries_df[['name', 'ref_area']]
 countries_df['name'] = countries_df.name.str.strip()
 countries_df.loc[len(countries_df)] = ['Antarctica', 14_200_000]  # from: https://en.wikipedia.org/wiki/Antarctica
 ```

 %% Cell type:markdown id:0004-5231db90b71a4c73fd2da36a3ca092775410c9eb4e21d28e7de691474ef tags:

 Hier kann ein EPSG-Code gewählt werden, der die Projektion definiert.
 Damit bestimmt sich wie groß die projizierte Fläche ist.

-%% Cell type:code id:0005-87f4bcf5159c2d37428cabbcd55f47bf4cd95bbb86ed15eced56f1b465b tags:
+%% Cell type:code id:0005-5b4adbf9e675eff3ed4415e452a13cfb859e618ba5637ea6a5290655f08 tags:

 ``` 
 countries = geo_countries.merge(countries_df, on='name')

 scaledown = False
 # epsg = 25832  # UTM Zone 32N, good for germany, does not work on global scale
 # epsg = 4087
 # epsg = 3395
 # epsg = 4326
 # ...

 # EPSG:3857 is the original projection CRS of the background map and thus the map corresponds to the computed projection error. So in this case we can scale geometries down to show true sizes of countries.
 epsg, scaledown = 3857, True

 countries = countries.to_crs(f'EPSG:{epsg}')

 # calculate relative difference of projected and reference areas
 countries['proj_area'] = countries.area / 1e6  # in km²
 ratio = countries['proj_area'] / countries['ref_area']
 countries['rel_diff_area'] = ((ratio - 1).round(2) * 100).astype(int)  # in percent

 if scaledown:
    # split MultiPolygons into separate Polygons, e. g. 1 MultiPolygon of Germany → 6 Polygons (6 rows with the same data)
-    countries = countries.explode()
+    countries = countries.explode(index_parts=False)

    # scale down each poly
    scale = 1 / np.sqrt(ratio)
    countries.geometry = [affinity.scale(g, s, s) for g, s in zip(countries.geometry, scale[countries.index])]

    # special handling for Antarctica, since it would be below the visible area
    max_antarctica_area = countries.loc[countries.name == 'Antarctica'].area.max()
    largest_antarctica_poly = countries.area == max_antarctica_area
    shifted_antarctica = affinity.translate(countries.loc[largest_antarctica_poly, 'geometry'].squeeze(), yoff=np.sqrt(max_antarctica_area)*7)
    countries.loc[largest_antarctica_poly, 'geometry'] = shifted_antarctica

    # reverse splitting, i. e. combine Polygons with the same name into a single MultiPolygon
    countries = countries.dissolve('name')

 # plot
 vmax = countries['rel_diff_area'].abs().quantile(0.97)
 countries.explore(
    column='rel_diff_area',
    vmin=-vmax, vmax=vmax, cmap='coolwarm',
    tooltip=['name', 'rel_diff_area', 'proj_area', 'ref_area'],
    popup=True,
    tiles="CartoDB voyager",
    attribution=countries.crs.name
 )
 ```

 %% Cell type:markdown id:0007-c166693366c20bec45d0a68d7bc44664690db85407f20dd66ceb12f7251 tags:

 Siehe auch:

 -   https://truesizeofcountries.com/
 -   https://www.thetruesize.com