Natura tiff check

[ekatef]

Closes #442 .

Continues #454 on the updated repo.

Changes proposed in this Pull Request

Discrepancy between the cutout and natura.tiff extents can lead to errors when generating renewable profiles. This PR is intended to prevent this by

• fixing natura.tiff generation (currently there seems to be a coordinates flip and CRS mismatch)
• adding a check to ensure that natura.tiff fully covers the cutout area and was generated using a proper CRS

Checklist

• I tested my contribution locally and it seems to work fine.
• Code and workflow changes are sufficiently documented.
• Newly introduced dependencies are added to envs/environment.yaml and envs/environment.docs.yaml.
• Changes in configuration options are added in all of config.default.yaml, config.tutorial.yaml, and test/config.test1.yaml.
• Changes in configuration options are also documented in doc/configtables/*.csv and line references are adjusted in doc/configuration.rst and doc/tutorial.rst.
• A note for the release notes doc/release_notes.rst is amended in the format of previous release notes, including reference to the requested PR.

[ekatef]

@davide-f, many thanks for checking!
I think, now I have answers on all the points related to natura.tiff generation.

I see the "natura_tiff_viz" is a python script; may you create is as a notebook, similarly to the others? I think this would be great from an usability perspective.

Totally agree! Such a visualisation tool could be very helpful to check cutout-natura.tiff correspondance. The notebook is there
https://github.com/ekatef/documentation/blob/geo_check/notebooks/code_dev_notebooks/check_natura_tiff.ipynb

Regarding the question you ose to Fabian, I believe keeping the lat-lon convenction coherent in the script should be the way to go
Great, then we should keep coordinates order for out_shape but remove a coordinate flip [::-1] when using out_shape.

Anyway, I'd expect the orange box to have an omothetic shape than the blue one; I think there is the need to fix that.
It would be nice if in the plot you could show the parallels/meridians (gridded map), so that it would be easier to understand.
My feeling is that (a) there may be a crs problem, (b) an inversion lat/lon, or (c) miscalculation of the boundaries.

The point is cutout and natura.tiff raster are being kept as rectangulars but were generated using different projections. Which means that we need to make a coordinates transform (on at least one of them) to be able overlay them on the same map. The grid is really helpful to highlight this effect: compare please the maps in the notebook.

Of course, I may be biased but for me it looks like now the problem is fixed as compared with an initial natura.tiff generation approach: for Africa or Central Asia

[davide-f]

Thanks @ekatef !

Regarding the projections that's unfortunate and it shall be fixed in future releases.
when you mean that natura and cutout are using different projections, what are the projections used for the two files?

For natura, there may have been an old refuse that may have led to using 3035 by mistake and that shall be fixed, if that's the case.

[ekatef]

@davide-f, thanks for the comments :)

Regarding the projections that's unfortunate and it shall be fixed in future releases.
when you mean that natura and cutout are using different projections, what are the projections used for the two files?

For natura, there may have been an old refuse that may have led to using 3035 by mistake and that shall be fixed, if that's the case.

Yes, CRS of the pre-built natura.tiff is EPSG:3035 while a freshly-generated version has ESRI:54009. The cutout derived from ERA5 nc-files implies EPSG:4326 (as geo_crs).

If I understand properly area_crs is needed to keep the correct area value free from distortions when merging shapes:

pypsa-earth/scripts/build_natura_raster.py

Line 150 in e35545c

unified_shape_file = unary_union(shape["geometry"])

...and generating the raster itself

pypsa-earth/scripts/build_natura_raster.py

Lines 187 to 188 in e35545c

	raster = ~geometry_mask(shapes.geometry, out_shape[::-1], transform)
	raster = raster.astype(rio.uint8)

Probably we could re-project the resulted raster and cut it to match exactly with to the cutout extend? However, I wouldn't say that raster approach to re-projection is the most elegant one (which seems to be connected with the raster format itself):

https://github.com/rasterio/rasterio/blob/main/examples/reproject.py

Have I understood your idea correctly?

[ekatef]

@davide-f, would be very grateful for your review.

I have left natura.tiff with area_crs for now. As far as I understand there should be no difference between the natura.tiff and cutout borders when we'll be using the global natura.tiff. For a "regional" natura.tiff version there will be parts of the raster uncovered with the cutout but I don't expect any issues due to this except a bit less effective memory usage. However, I think that's better than adding re-projection manipulations which can be a bit tricky to deal with. Do you agree?

[davide-f]

@davide-f, would be very grateful for your review.

I have left natura.tiff with area_crs for now. As far as I understand there should be no difference between the natura.tiff and cutout borders when we'll be using the global natura.tiff. For a "regional" natura.tiff version there will be parts of the raster uncovered with the cutout but I don't expect any issues due to this except a bit less effective memory usage. However, I think that's better than adding re-projection manipulations which can be a bit tricky to deal with. Do you agree?

Perfect! I completely agree.
As long as the calculations are right, better not to overmanipulate things.

@ekatef do you think this PR is ready?

[ekatef]

@davide-f Thank you :)
Yes, I think it's ready for the review

[davide-f]

@ekatef I created the new world natura.tiff, may you check it with your scripts and compare to a cutout?
Link to the new natura file (from PR 470). https://drive.google.com/file/d/1Owl38VQj1_pgwcGFvfxpX94moxtpcqSR/view?usp=sharing

[ekatef]

@davide-f, fantastic news :) Thank you so much!
Have loaded the global natura.tiff. Will be happy to play with it and notify you about the results as soon as I can.

[ekatef]

@davide-f, unfortunately it looks like there are some troubles with applying the global natura.tiff. Coordinates transformation to any cylindrical projection ("EPSG:4326" or "EPSG:32663") gives the outer natura.tiff boundary as POLYGON ((Infinity Infinity, Infinity Infinity, Infinity Infinity, Infinity Infinity, Infinity Infinity)) while building renewable profiles with Central Asia cutout leads to an error Point outside of projection domain (the full listing attached bellow).

That's a bit strange as the original raster boundary is POLYGON ((-359000 -27057000, -359000 9020100, -18040100 9020100, -18040100 -27057000, -359000 -27057000)) which gives 36077100 to 17681100 with an aspect rasio 2/1 how it should be...

Just to be sure about coordinates definition:

https://github.com/davide-f/pypsa-earth/blob/036bcc8d67ed423743daf7b0549887b382e7b526/scripts/build_natura_raster.py#L98-L101

Should be here the limiting range be +/-180 deg the same both for x and y coordinates?

# ----------------------------------------------------
# error listing
# ----------------------------------------------------
INFO:snakemake.logging:rule build_renewable_profiles:
    input: networks/base.nc, resources/natura.tiff, data/copernicus/PROBAV_LC100_global_v3.0.1_2019-nrt_Discrete-Classification-map_EPSG-4326.tif, data/gebco/GEBCO_2021_TID.nc, resources/shapes/country_shapes.geojson, resources/shapes/offshore_shapes.geojson, data/hydro_capacities.csv, data/eia_hydro_annual_generation.csv, resources/powerplants.csv, resources/bus_regions/regions_onshore.geojson, cutouts/casr-2018-era5.nc
    output: resources/renewable_profiles/profile_onwind.nc
    log: logs/build_renewable_profile_onwind.log
    jobid: 8
    benchmark: benchmarks/build_renewable_profiles_onwind
    reason: Missing output files: resources/renewable_profiles/profile_onwind.nc; Input files updated by another job: networks/base.nc, resources/bus_regions/regions_onshore.geojson, resources/shapes/offshore_shapes.geojson, resources/shapes/country_shapes.geojson, resources/powerplants.csv
    wildcards: technology=onwind
    resources: tmpdir=/var/folders/qn/vpndfm21795ckkq89np1ckp40000gn/T, mem_mb=20000

INFO:snakemake.logging:
Traceback (most recent call last):
  File "~Documents/_github_/pypsa-earth/.snakemake/scripts/tmpsjdn5q29.build_renewable_profiles.py", line 425, in <module>
    natura_gejson = rio.warp.transform_geom(
  File "~opt/miniconda3/envs/pypsa-africa/lib/python3.9/site-packages/rasterio/env.py", line 387, in wrapper
    return f(*args, **kwds)
  File "~opt/miniconda3/envs/pypsa-africa/lib/python3.9/site-packages/rasterio/env.py", line 612, in wrapper
    return f(*args, **kwds)
  File "~opt/miniconda3/envs/pypsa-africa/lib/python3.9/site-packages/rasterio/warp.py", line 109, in transform_geom
    return _transform_geom(
  File "rasterio/_warp.pyx", line 125, in rasterio._warp._transform_geom
  File "rasterio/_warp.pyx", line 76, in rasterio._warp._transform_single_geom
  File "rasterio/_err.pyx", line 215, in rasterio._err.exc_wrap_pointer
rasterio._err.CPLE_AppDefinedError: Point outside of projection domain
[Thu Oct 13 22:23:50 2022]
INFO:snakemake.logging:[Thu Oct 13 22:23:50 2022]
Error in rule build_renewable_profiles:
    jobid: 8
    output: resources/renewable_profiles/profile_onwind.nc
    log: logs/build_renewable_profile_onwind.log (check log file(s) for error message)

ERROR:snakemake.logging:Error in rule build_renewable_profiles:
    jobid: 8
    output: resources/renewable_profiles/profile_onwind.nc
    log: logs/build_renewable_profile_onwind.log (check log file(s) for error message)

RuleException:
CalledProcessErrorin line 341 of ~Documents/_github_/pypsa-earth/Snakefile:
Command 'set -euo pipefail;  ~opt/miniconda3/envs/pypsa-africa/bin/python3.9 ~Documents/_github_/pypsa-earth/.snakemake/scripts/tmpsjdn5q29.build_renewable_profiles.py' returned non-zero exit status 1.
  File "~Documents/_github_/pypsa-earth/Snakefile", line 341, in __rule_build_renewable_profiles
  File "~opt/miniconda3/envs/pypsa-africa/lib/python3.9/concurrent/futures/thread.py", line 58, in run
ERROR:snakemake.logging:RuleException:
CalledProcessErrorin line 341 of ~Documents/_github_/pypsa-earth/Snakefile:
Command 'set -euo pipefail;  ~opt/miniconda3/envs/pypsa-africa/bin/python3.9 ~Documents/_github_/pypsa-earth/.snakemake/scripts/tmpsjdn5q29.build_renewable_profiles.py' returned non-zero exit status 1.
  File "~Documents/_github_/pypsa-earth/Snakefile", line 341, in __rule_build_renewable_profiles
  File "~opt/miniconda3/envs/pypsa-africa/lib/python3.9/concurrent/futures/thread.py", line 58, in run
Shutting down, this might take some time.
WARNING:snakemake.logging:Shutting down, this might take some time.
Exiting because a job execution failed. Look above for error message
ERROR:snakemake.logging:Exiting because a job execution failed. Look above for error message
Complete log: .snakemake/log/2022-10-13T220051.621624.snakemake.log
WARNING:snakemake.logging:Complete log: .snakemake/log/2022-10-13T220051.621624.snakemake.log

[davide-f]

Thanks katia! I'll investigate the issue; I reposted the comment in the other PR.
Let's move the discussion there

[davide-f]

I'm wondering: rather than testing the cutout to be contained into natura, shouldn't we check:

• the "world_shape" to be contained into natura
• the "world_shape" to be contained into cutout

Having cutout contained into natura is not necessarily what we need.
What do you think?

[ekatef]

I'm wondering: rather than testing the cutout to be contained into natura, shouldn't we check:

• the "world_shape" to be contained into natura
• the "world_shape" to be contained into cutout

Having cutout contained into natura is not necessarily what we need. What do you think?

Actually, this PR was intended to fix a mismatch between natura.tiff and cutout areas. Such a mismatch may be difficult to detect but if being present it leads to issues like reported here #411 (comment).

Agree that it may be worth to re-think checking of the GIS inputs. My feeling is that the country-cutout correspondence check is definitely needed.

As for natura.tiff, the idea was to ensure that the raster file was generated correctly. Since build_natura_raster utilises only the cutout and does not include the country shape, in case of mismatch between the natura raster and the cutout, we may suspect that something could went wrong.

From the perspective of the model interfaces the logic your suggest may be more natural. E.g. we supply the intersecting cutout and natura.tiff while both are fine for a selected country, and the workflow should be happy with that. But the proposed check will not have any objections in that case just issuing a warning which indicates that the situation is a bit unusual.

[ekatef]

By the way, the problem with renewable_profiles_generation for China was linked to data loss in the raw ERA5 data. Many thanks @davide-f for sharing experiences in this regard :-)

My feeling is that probably it may be worth to tackle both cutout issues (location and data validity) in a separate PR:

• check that the cutout corresponds to the requested region;
• check is there are any data loss in the cutout data;
• fix NaN when processing the capacity_factor array.

@davide-f, what do you think?

[ekatef]

@pz-max, could you please check if merging of this PR went properly?
I may be wrong but I'm not able to see in the main any changes I have introduced here...

[pz-max]

@pz-max, could you please check if merging of this PR went properly? I may be wrong but I'm not able to see in the main any changes I have introduced here...

This is super weird... the main doesn't show it. Once, I click on "show difference to the previous version, it will show up in the original & the previous version...

[ekatef]

This is super weird... the main doesn't show it. Once, I click on "show difference to the previous version, it will show up in the original & the previous version...

@pz-max, thank you for checking! I have done a local copy of the PR branch. In case it'll not show-up in some days, I'll just quickly create a new version.

Software is generally unpredictable... (c) :D