From 8538233bba5d0c24f86e8c296198546657ca5d6b Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Thu, 7 May 2020 14:45:14 +0200 Subject: [PATCH] Fix bug due to invalid cluster bus region geometries These geometries are apparently invalid due to self-crossing or self-touching polygons. The geometries are created by pypsa-eur/scripts/cluster_network.py but appear to be valid before being written to file. They are only valid after being read back in from file. This seems to indicate some numerical issue relating to file reading and writing. Now the geometries are cleaned after being read in. --- scripts/build_clustered_population_layouts.py | 4 ++-- scripts/build_heat_demand.py | 4 +++- scripts/build_solar_thermal_profiles.py | 9 +++------ scripts/build_temperature_profiles.py | 4 +++- scripts/helper.py | 15 +++++++++++++++ 5 files changed, 26 insertions(+), 10 deletions(-) create mode 100644 scripts/helper.py diff --git a/scripts/build_clustered_population_layouts.py b/scripts/build_clustered_population_layouts.py index e7f4db5a..94cbd133 100644 --- a/scripts/build_clustered_population_layouts.py +++ b/scripts/build_clustered_population_layouts.py @@ -3,7 +3,7 @@ import geopandas as gpd import xarray as xr import pandas as pd import atlite - +import helper cutout = atlite.Cutout(snakemake.config['atlite']['cutout_name'], @@ -14,7 +14,7 @@ clustered_busregions_as_geopd = gpd.read_file(snakemake.input.regions_onshore).s clustered_busregions = pd.Series(clustered_busregions_as_geopd.geometry, index=clustered_busregions_as_geopd.index) - +helper.clean_invalid_geometries(clustered_busregions) I = cutout.indicatormatrix(clustered_busregions) diff --git a/scripts/build_heat_demand.py b/scripts/build_heat_demand.py index 1ed57ea4..865f12bd 100644 --- a/scripts/build_heat_demand.py +++ b/scripts/build_heat_demand.py @@ -4,7 +4,7 @@ import atlite import pandas as pd import xarray as xr import scipy as sp - +import helper if 'snakemake' not in globals(): from vresutils import Dict @@ -26,6 +26,8 @@ clustered_busregions_as_geopd = gpd.read_file(snakemake.input.regions_onshore).s clustered_busregions = pd.Series(clustered_busregions_as_geopd.geometry, index=clustered_busregions_as_geopd.index) +helper.clean_invalid_geometries(clustered_busregions) + I = cutout.indicatormatrix(clustered_busregions) diff --git a/scripts/build_solar_thermal_profiles.py b/scripts/build_solar_thermal_profiles.py index a3a171d9..c26266aa 100644 --- a/scripts/build_solar_thermal_profiles.py +++ b/scripts/build_solar_thermal_profiles.py @@ -1,15 +1,10 @@ -import sys - -#override atlite -sys.path = ["/home/vres/data/tom/lib/atlite"] + sys.path - import geopandas as gpd import atlite import pandas as pd import xarray as xr import scipy as sp - +import helper if 'snakemake' not in globals(): from vresutils import Dict @@ -33,6 +28,8 @@ clustered_busregions_as_geopd = gpd.read_file(snakemake.input.regions_onshore).s clustered_busregions = pd.Series(clustered_busregions_as_geopd.geometry, index=clustered_busregions_as_geopd.index) +helper.clean_invalid_geometries(clustered_busregions) + I = cutout.indicatormatrix(clustered_busregions) diff --git a/scripts/build_temperature_profiles.py b/scripts/build_temperature_profiles.py index 28a4b858..a55bd606 100644 --- a/scripts/build_temperature_profiles.py +++ b/scripts/build_temperature_profiles.py @@ -4,7 +4,7 @@ import atlite import pandas as pd import xarray as xr import scipy as sp - +import helper if 'snakemake' not in globals(): from vresutils import Dict @@ -27,6 +27,8 @@ clustered_busregions_as_geopd = gpd.read_file(snakemake.input.regions_onshore).s clustered_busregions = pd.Series(clustered_busregions_as_geopd.geometry, index=clustered_busregions_as_geopd.index) +helper.clean_invalid_geometries(clustered_busregions) + I = cutout.indicatormatrix(clustered_busregions) diff --git a/scripts/helper.py b/scripts/helper.py new file mode 100644 index 00000000..6f060db4 --- /dev/null +++ b/scripts/helper.py @@ -0,0 +1,15 @@ + +import logging +logger = logging.getLogger(__name__) + +#https://stackoverflow.com/questions/20833344/fix-invalid-polygon-in-shapely +#https://stackoverflow.com/questions/13062334/polygon-intersection-error-in-shapely-shapely-geos-topologicalerror-the-opera +#https://shapely.readthedocs.io/en/latest/manual.html#object.buffer +def clean_invalid_geometries(geometries): + """Fix self-touching or self-crossing polygons; these seem to appear +due to numerical problems from writing and reading, since the geometries +are valid before being written in pypsa-eur/scripts/cluster_network.py""" + for i,p in geometries.items(): + if not p.is_valid: + logger.warning(f'Clustered region {i} had an invalid geometry, fixing using zero buffer.') + geometries[i] = p.buffer(0)