Harmonize CRS in the workflow

config.default.yaml

@@ -24,7 +24,7 @@ summary_dir: results
 snapshots:
   start: "2013-01-01"
   end: "2014-01-01"
-  closed: "left" # end is not inclusive
+  inclusive: "left" # end is not inclusive
 
 enable:
   retrieve_databundle: true  #  Recommended 'true', for the first run. Otherwise data might be missing.
@@ -35,6 +35,12 @@ enable:
   # requires cds API key https://cds.climate.copernicus.eu/api-how-to
   # More information https://atlite.readthedocs.io/en/latest/introduction.html#datasets
 
+# definition of the Coordinate Reference Systems
+crs:
+  geo_crs: EPSG:4326  # general geographic projection, not used for metric measures. "EPSG:4326" is the standard used by OSM and google maps
+  distance_crs: EPSG:3857  # projection for distance measurements only. Possible recommended values are "EPSG:3857" (used by OSM and Google Maps)
+  area_crs: ESRI:54009  # projection for area measurements only. Possible recommended values are Global Mollweide "ESRI:54009"
+
 # download_osm_data_nprocesses: 10  # (optional) number of threads used to download osm data
 
 augmented_line_connection:

config.tutorial.yaml

@@ -80,7 +80,7 @@ summary_dir: results
 snapshots:
   start: "2013-03-1"
   end: "2013-03-14"
-  closed: "left" # end is not inclusive
+  inclusive: "left" # end is not inclusive
 
 enable:
   # prepare_links_p_nom: false
@@ -91,6 +91,12 @@ enable:
   build_cutout: false
   build_natura_raster: true  # If True, then build_natura_raster can be run
 
+# definition of the Coordinate Reference Systems
+crs:
+  geo_crs: EPSG:4326  # general geographic projection, not used for metric measures. "EPSG:4326" is the standard used by OSM and google maps
+  distance_crs: EPSG:3857  # projection for distance measurements only. Possible recommended values are "EPSG:3857" (used by OSM and Google Maps)
+  area_crs: ESRI:54009  # projection for area measurements only. Possible recommended values are Global Mollweide "ESRI:54009"
+
 electricity:
   voltages: [220., 300., 380.]
   co2limit: 7.75e+7 # 0.05 * 3.1e9*0.5

doc/release_notes.rst

@@ -45,6 +45,8 @@ Upcoming Release
 
 * Add latitude_optimal to get optimal solar orientation by default `Commit 1b2466b <https://github.com/pypsa-meets-africa/pypsa-africa/commit/de7d32be8807e4fc42486a60184f45680612fd46>`_
 
+* Harmonize CRSs by options `PR #356 <https://github.com/pypsa-meets-africa/pypsa-africa/pull/356>`_
+
 
 PyPSA-Africa 0.0.2 (6th April 2022)
 =====================================

scripts/base_network.py

@@ -214,20 +214,18 @@ def _set_countries_and_substations(n):
     buses = n.buses
 
     countries = snakemake.config["countries"]
-    country_shapes = (
-        gpd.read_file(snakemake.input.country_shapes)
-        .set_index("name")["geometry"]
-        .set_crs(4326)
-    )
+    country_shapes = gpd.read_file(snakemake.input.country_shapes).set_index("name")[
+        "geometry"
+    ]
     offshore_shapes = unary_union(
-        gpd.read_file(snakemake.input.offshore_shapes)["geometry"].set_crs(4326)
+        gpd.read_file(snakemake.input.offshore_shapes)["geometry"]
     )
 
     bus_locations = buses
     bus_locations = gpd.GeoDataFrame(
         bus_locations,
         geometry=gpd.points_from_xy(bus_locations.x, bus_locations.y),
-        crs=4326,
+        crs=country_shapes.crs,  # the workflow sets the the same crs for buses and shapes
     )
     # Check if bus is in shape
     offshore_b = bus_locations.within(offshore_shapes)

scripts/build_natura_raster.py

@@ -62,6 +62,8 @@
 _logger = logging.getLogger(__name__)
 _logger.setLevel(logging.INFO)
 
+CUTOUT_CRS = "EPSG:4326"
+
 
 def get_fileshapes(list_paths, accepted_formats=(".shp",)):
     "Function to parse the list of paths to include shapes included in folders, if any"
@@ -88,7 +90,7 @@ def determine_cutout_xXyY(cutout_name, out_logging):
     if out_logging:
         _logger.info("Stage 1/5: Determine cutout boundaries")
     cutout = atlite.Cutout(cutout_name)
-    assert cutout.crs.to_epsg() == 4326
+    assert cutout.crs == CUTOUT_CRS
     x, X, y, Y = cutout.extent
     dx, dy = cutout.dx, cutout.dy
     return [x - dx / 2.0, X + dx / 2.0, y - dy / 2.0, Y + dy / 2.0]
@@ -104,7 +106,7 @@ def get_transform_and_shape(bounds, res, out_logging):
     return transform, shape
 
 
-def unify_protected_shape_areas(inputs, out_logging):
+def unify_protected_shape_areas(inputs, area_crs, out_logging):
     """
     Iterates thorugh all snakemake rule inputs and unifies shapefiles (.shp) only.
 
@@ -113,7 +115,7 @@ def unify_protected_shape_areas(inputs, out_logging):
 
     Returns
     -------
-    unified_shape : GeoDataFrame with a unified "multishape" (crs=3035)
+    unified_shape : GeoDataFrame with a unified "multishape"
 
     """
     import pandas as pd
@@ -133,16 +135,14 @@ def unify_protected_shape_areas(inputs, out_logging):
     for i in shp_files:
         shape = gpd.GeoDataFrame(
             pd.concat([gpd.read_file(i) for i in shp_files])
-        ).to_crs(3035)
+        ).to_crs(area_crs)
 
     # Removes shapely geometry with null values. Returns geoseries.
     shape = shape["geometry"][shape["geometry"].is_valid]
 
-    # Create Geodataframe with crs(3035)
-    shape = gpd.GeoDataFrame(shape)
-    shape = shape.rename(columns={0: "geometry"}).set_geometry(
-        "geometry"
-    )  # .set_crs(3035)
+    # Create Geodataframe with crs
+    shape = gpd.GeoDataFrame(shape, crs=area_crs)
+    shape = shape.rename(columns={0: "geometry"}).set_geometry("geometry")
 
     # Unary_union makes out of i.e. 1000 shapes -> 1 unified shape
     if out_logging:
@@ -152,7 +152,7 @@ def unify_protected_shape_areas(inputs, out_logging):
         _logger.info(
             "Stage 3/5: Unify protected shape area. Step 3: Set geometry of unified shape"
         )
-    unified_shape = gpd.GeoDataFrame(geometry=[unified_shape_file]).set_crs(3035)
+    unified_shape = gpd.GeoDataFrame(geometry=[unified_shape_file], crs=area_crs)
 
     return unified_shape
 
@@ -165,19 +165,22 @@ def unify_protected_shape_areas(inputs, out_logging):
         snakemake = mock_snakemake("build_natura_raster")
     configure_logging(snakemake)
 
+    # get crs
+    area_crs = snakemake.config["crs"]["area_crs"]
+
     out_logging = True
     inputs = snakemake.input
     cutouts = inputs.cutouts
     shapefiles = get_fileshapes(inputs)
     xs, Xs, ys, Ys = zip(
         *(determine_cutout_xXyY(cutout, out_logging=out_logging) for cutout in cutouts)
     )
-    bounds = transform_bounds(4326, 3035, min(xs), min(ys), max(Xs), max(Ys))
+    bounds = transform_bounds(CUTOUT_CRS, area_crs, min(xs), min(ys), max(Xs), max(Ys))
     transform, out_shape = get_transform_and_shape(
         bounds, res=100, out_logging=out_logging
     )
     # adjusted boundaries
-    shapes = unify_protected_shape_areas(shapefiles, out_logging=out_logging)
+    shapes = unify_protected_shape_areas(shapefiles, area_crs, out_logging=out_logging)
 
     if out_logging:
         _logger.info("Stage 4/5: Mask geometry")
@@ -193,7 +196,7 @@ def unify_protected_shape_areas(inputs, out_logging):
         dtype=rio.uint8,
         count=1,
         transform=transform,
-        crs=3035,
+        crs=area_crs,
         compress="lzw",
         width=raster.shape[1],
         height=raster.shape[0],