Adding Iceland to analysis also adds random islands

[brynpickering]

This issue relates to adding Iceland for use in a downstream dependency, Euro-Calliope (see https://github.com/timtroendle/euro-calliope/issues/15).

Including Iceland here is relatively straightforward, it just requires:

1. Updating the snakemake config to include Iceland in the list of countries (inc. at every resolution), and to extend the spatial bounds to cover Iceland.
2. Updating the renewables.ninja capacity factor datasets, now that there are more points to simulate.

There's only one issue: on creating the ninja input configs I noticed that the extended spatial bounds of the problem leads to a bunch of island territories ending up in the model:

These same islands will end up elsewhere in the analysis, since the spatial bounds are used to set the study area in e.g. generating shapefiles. Thoughts on the easiest way to exclude these @timtroendle, or do you imagine they won't be a problem?

[timtroendle]

I think this is a problem as the demand data does not include these islands (entso-e, 2016). However, as this is the case, maybe it would be even good to have a special case treatment here that removes Madeira and the Azores explicitly, not implicitly like before.

entso-e. (2016). Specific national considerations. https://www.entsoe.eu/Documents/Publications/Statistics/Specific_national_considerations.pdf

[brynpickering]

I wonder if we could provide the bounds of the system as a shapefile, instead of bounds of a rectangle? Or does that over-complicate the cutting down of input data (which currently uses those bounds)?

[timtroendle]

I think the code would not change much if we were to use a shapefile (basically a complex polygon I guess, no?) rather than a box. The configuration would become more complicated though. Instead of 4 values you would need a long textfile describing the polygon (assuming it is stored as a GeoJSON rather than as a shapefile).

I think I'd rather take the current box, and then stamp our Madeira and the Azores. The coordinates of the islands can be hard coded. The stamping could be controlled by one or two more flags in the config:

scope:
    countries:
        [...]
    bounds:
        x_min: [...]
        [...]
        y_max: [...]
        include-madeira: True
        include-azores: True

Those are weirdly specific configuration parameters, but it looks like a pragmatic solution to me.

[brynpickering]

How about being able to define exclusion zones?

scope:
    countries:
        [...]
    bounds:
        x_min: [...]
        [...]
        y_max: [...]
    exclusion_zones:
        atlantic_islands: 
            x_min: [...]
            [...]
            y_max: [...]

It could then be any arbitrary number of exclusion zones that are named by the user for convenience (atlantic_zones could just as easily be random_stuff for all the workflow cares)

[timtroendle]

Much better!

[brynpickering]

OK, there are a few locations were bounds are used in the analysis:

1. in the clipping of rasters with rio clip --bounds {x_min},{x_max},{y_min},{y_max} {input} {output} CLI

This might need to be updated with clipping relative to a raster which has the 'holes' of the exclusion zones already in it: rio clip --like {some_raster.tif} {input} {output}. I'm unsure how to implement this. Does rasterio respect 'holes' in a raster that it's clipping data against, or does it just see the outline of the raster, and nothing about the NoValue pixels within it? If it's not possible, would it instead be possible to nullify data in the exclusion zones after clipping?

2. in the getting of nuts/gadm data with shapely.geometry.box(...) to get a study area.

Seems quite easy to wrap this in a Polygon shape, holes and all: shapely.geometry.Polygon(shapely.geometry.box({bounds}), holes=[shapely.geometry.box({i}) for i in exclusion_bounds])

[timtroendle]

1. I don't think it's necessary to change this. There is no need for the rasters to have holes. The important bit is that the shapes are correct, because they are used to "stamp" rasters. That's why...

2. ... is important.

3. There is, however, another part of the workflow that may need to be touched when changing the bounds: the bounds of the elevation data. It would be better if those bounds were derived from the configuration values instead of being hard coded. This has nothing to do with the "holes" discussed in this issue though -- it just came into my mind.

[brynpickering]

Great, well 2. is easy - I'll create a PR on that soon. As to SRTM tiles, is it possible to find the spatial bounds of tiles by their numbers? On another note, is there more hardcoding going on here with the 59.5 and 60 values? https://github.com/timtroendle/possibility-for-electricity-autarky/blob/c7bf64c68ab1810c42f9b42a5624ad2419920891/rules/data-preprocessing.smk#L286-L287

[timtroendle]

No, the 59.5 and 60 values are fine. That's just because SRTM is only available up until latitude 60°.

For the SRTM files it's likely possible to map bounds to their numbers. I simply did not bother back then.

[brynpickering]

SRTM: srtm_X_Y.zip where X = (lon + 180)/5 + 1 and Y = (55 - lat)/5 + 1. lon and lat are the bottom left corner of a 5*5 degree grid.