Natura tiff check

[ekatef]

Closes #442 .

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 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]

An initial idea of the proposed fix implied that natura.tiff raster should cover the whole cutout area. However, visual checks have resulted in the following pictures (blue for cutout, orange for natura.tiff):

Africa

Central Asia

It appears that boundaries of the natura.tiff raster may be not perfectly correspond to the cutout extend. I suspect that x and y coordinates could be swapped when defining the shape variable which is used by raster generation and contains number of rows and columns.

First, the shape is defined like that
https://github.com/pypsa-meets-africa/pypsa-africa/blob/1065e4641b0fd1519fe4fcef999d2522cfcc449c/scripts/build_natura_raster.py#L104

Then it's reversed when using it when specifying raster
https://github.com/pypsa-meets-africa/pypsa-africa/blob/1065e4641b0fd1519fe4fcef999d2522cfcc449c/scripts/build_natura_raster.py#L187

Changing the order in the shape definition seems to fix the issue
https://github.com/ekatef/pypsa-africa/blob/6c507c60ba14535cc6e8c44cc2aede71fa6c04ad/scripts/build_natura_raster.py#L107-L108

The cutout-natura.tiff overlay after applying this change looks like that:

However, I'm not quite sure that have completely understood the logic behind the raster transformation.

@davide-f, @pz-max, would be very grateful if you could please have a look on this.

[pz-max]

@ekatef thanks for this amazing report 🥇 Yes, a check of all of the blue area is covered by orange would be great. Ideally, the nature raster should match the cutout which seemingly did not work. You fix improved things but it seems there is some unnecessary stuff happening? I think the aim is to have a little buffer/ more of the natura area compared to the cutout (=orange being slightly larger than blue).

Regarding the raster transformation, I am also not quite sure maybe @euronion has an idea

[ekatef]

@ekatef thanks for this amazing report 🥇 Yes, a check of all of the blue area is covered by orange would be great. Ideally, the nature raster should match the cutout which seemingly did not work. You fix improved things but it seems there is some unnecessary stuff happening? I think the aim is to have a little buffer/ more of the natura area compared to the cutout (=orange being slightly larger than blue).

Regarding the raster transformation, I am also not quite sure maybe @euronion has an idea

Many thanks, @pz-max!

If I understand properly, a safety margin is guaranteed by theout_shape definition:

https://github.com/pypsa-meets-africa/pypsa-africa/blob/1065e4641b0fd1519fe4fcef999d2522cfcc449c/scripts/build_natura_raster.py#L102-L103

So, both the raster dimensions are on one res higher as compared to the cutout bounds, right?

What I don't understand is why the left upper corner of the cutout and the natura-raster do not match each other for Africa. According to raster generation in build_natura_raster it looks like they should...

[pz-max]

Yes, it raises question marks. Let's discuss that tomorrow at the weekly. Super interesting and important issue anyways. Thanks for spotting this 🥇

[ekatef]

Update: it looks like that the reason of the shape discrepancy between the natura.tiff raster and the cutout for Africa may by a difference between CRSs of the natura.tiff and the area_crs.

It seemsnatura.tiff available via the bundle_config may have been generated using EPSG:3035 (Lambert azimuthal) instead of ESRI:54009 (Mollweide pseudocylindrical). Actually, that is this CRS discrepancy not coverage issues which lead to the error appeared in #411.

[ekatef]

Have checked natura.tiff generation for Africa, Central Asia and Argentina (to have a test for Southern hemisphere) and it seems to work fine.

Africa

Argentina

It could also make sense to update the natura.tiff provided with the model.

[ekatef]

I think that's ready for the review.
@pz-max, @davide-f, would be very grateful if you could please look into it.

[euronion]

@ekatef thanks for this amazing report 🥇 Yes, a check of all of the blue area is covered by orange would be great. Ideally, the nature raster should match the cutout which seemingly did not work. You fix improved things but it seems there is some unnecessary stuff happening? I think the aim is to have a little buffer/ more of the natura area compared to the cutout (=orange being slightly larger than blue).

Regarding the raster transformation, I am also not quite sure maybe @euronion has an idea

Better to direct questions regarding the transform to @FabianHofmann .

In general note that with newer atlite versions you don't need to create the raster if you don't want to, but can directly use the WDPA data :)

[FabianHofmann]

@ekatef are you sure about switching the order of x y? Normally in GIS applications, the y axis is the first dimension (as in a matrix) and x the second which is indeed confusing.

[ekatef]

Better to direct questions regarding the transform to @FabianHofmann .

In general note that with newer atlite versions you don't need to create the raster if you don't want to, but can directly use the WDPA data :)

@euronion, many thanks for looking into it! Good to know about changes in the atlite approach. It sounds like the new one will be more convenient :)

[ekatef]

@FabianHofmann, many thanks for the review.

I have been trying to keep the logic of the original code where a coordinates swap is applied to the out_shape when using it for raster generation:

https://github.com/pypsa-meets-africa/pypsa-africa/blob/1065e4641b0fd1519fe4fcef999d2522cfcc449c/scripts/build_natura_raster.py#L187

If I understand properly, the intention behind this approach was to have a regular cartesian coordinates order (x, y) in out_shape itself and swap it just before calculations where a reverse order is needed (n_rows = y, n_columns = x). But probably it would be better to keep a "geophysical" coordinates order in out_shape (y = lat, x = long) and remove [::-1]-swap instead. What do you think?

Generally, absolutely agree that it's always worth to be aware of "x for longitude, y for latitude" when working with coordinates for GIS-related problems...

[ekatef]

The testing procedure for an arbitrary region may be reproduces as follows:

1. Set config for a region/country of interest
2. Build cutout (mind Copernicus API...)
3. Build natura.tiff
4. Run a visualizer (./notebooks/natura_tiff_viz.py) [not sure if it makes sense to merge the visualising script as well]

@pz-max, will be very grateful for your review :)

[davide-f]

The testing procedure for an arbitrary region may be reproduces as follows:

1. Set config for a region/country of interest
2. Build cutout (mind Copernicus API...)
3. Build natura.tiff
4. Run a visualizer (./notebooks/natura_tiff_viz.py) [not sure if it makes sense to merge the visualising script as well]

@pz-max, will be very grateful for your review :)

Awesome @ekatef !
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.
I totally agree on the procedure that you suggested.

Regarding the question you ose to Fabian, I believe keeping the lat-lon convenction coherent in the script should be the way to go
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.
What is the status of the code with these latest changes?

[pz-max]

@ekatef , can we solve this today? To my understanding from the weekly meeting, I could just review this PR as it appears to work already.

[ekatef]

@davide-f, thank you so much for looking into it! I'll check all the points you have raised.

@pz-max, not sure I'll be able to finalize it today. I'd expect that some additional work with data will be needed to double-check the proposed approach.

But it'll be not difficult even to re-create this PR if needed as the most work here is work with data not code. So, if you are intended to rebrand today I think it absolutely makes sense to proceed with it :)

[pz-max]

@ekatef thanks for the feedback. Keep your local changes to be contributed from a fresh clone in future.
I am closing this PR now 💯

[pz-max]

@ekatef , rebranding is done. Feel free to work on the PR from a fresh clone or even fresh fork.

[ekatef]

@pz-max, many thanks for the great news :) Amazing work!

Regarding the new fork: it looks like I need to remove the previous fork to create a new one, right?

[pz-max]

If you do a clone from the fork.. How large is the repo?
A fresh fork is about 3.8MB

[ekatef]

A fresh clone from a fresh fork in my case is 3.1MB.

[pz-max]

Even better :)