search

PyPSA / atlite

Atlite: A Lightweight Python Package for Calculating Renewable Power Potentials and Time Series
https://atlite.readthedocs.io
254 stars 86 forks source link

Atlite errors with ESRI:540060 reprojections and Fiji #301

Open davide-f opened 1 year ago

davide-f commented 1 year ago

Description

When using Atlite on regions across borders of [+-180 lon /+-90° lan] GDAL returns error rasterio._err.CPLE_NotSupportedError: Cannot find coordinate operations from 'ESRI:54009' to 'EPSG:4326'.

The error is reproducible executing the PyPSA-Earth workflow on "FJ".

Error Message

Had to go 1 folder(s) up.
INFO:__main__:correction_factor is set as 0.8855
INFO:__main__:Calculate landuse availabilities...
INFO:rasterio._env:GDAL signalled an error: err_no=1, msg='Too many points (529 out of 529) failed to transform, unable to compute output bounds.'
WARNING:rasterio._env:CPLE_AppDefined in Unable to compute source region for output window 2479,5616,2479,1873, skipping.
INFO:rasterio._env:GDAL signalled an error: err_no=1, msg='Too many points (529 out of 529) failed to transform, unable to compute output bounds.'
WARNING:rasterio._env:CPLE_AppDefined in Unable to compute source region for output window 2479,5616,2479,1873, skipping.
INFO:rasterio._env:GDAL signalled an error: err_no=1, msg='Too many points (529 out of 529) failed to transform, unable to compute output bounds.'
WARNING:rasterio._env:CPLE_AppDefined in Unable to compute source region for output window 0,6332,2730,3166, skipping.
INFO:rasterio._env:GDAL signalled an error: err_no=1, msg='Too many points (529 out of 529) failed to transform, unable to compute output bounds.'
WARNING:rasterio._env:CPLE_AppDefined in Unable to compute source region for output window 0,9498,2730,3167, skipping.
INFO:rasterio._env:GDAL signalled an error: err_no=1, msg='Too many points (529 out of 529) failed to transform, unable to compute output bounds.'
WARNING:rasterio._env:CPLE_AppDefined in Unable to compute source region for output window 2730,11081,2730,1584, skipping.
INFO:rasterio._env:GDAL signalled an error: err_no=1, msg='Too many points (529 out of 529) failed to transform, unable to compute output bounds.'
WARNING:rasterio._env:CPLE_AppDefined in Unable to compute source region for output window 0,6332,2730,3166, skipping.
INFO:rasterio._env:GDAL signalled an error: err_no=1, msg='Too many points (529 out of 529) failed to transform, unable to compute output bounds.'
WARNING:rasterio._env:CPLE_AppDefined in Unable to compute source region for output window 0,9498,2730,3167, skipping.
INFO:rasterio._env:GDAL signalled an error: err_no=1, msg='Too many points (529 out of 529) failed to transform, unable to compute output bounds.'
WARNING:rasterio._env:CPLE_AppDefined in Unable to compute source region for output window 2730,11081,2730,1584, skipping.
INFO:rasterio._env:GDAL signalled an error: err_no=6, msg="Cannot find coordinate operations from `ESRI:54009' to `EPSG:4326'"
multiprocessing.pool.RemoteTraceback: 
"""
Traceback (most recent call last):
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.8/multiprocessing/pool.py", line 125, in worker
    result = (True, func(*args, **kwds))
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.8/multiprocessing/pool.py", line 48, in mapstar
    return list(map(*args))
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.8/site-packages/atlite/gis.py", line 673, in _process_func
    return shape_availability_reprojected(shapes.loc[[i]], *args)[0]
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.8/site-packages/atlite/gis.py", line 377, in shape_availability_reprojected
    return rio.warp.reproject(
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.8/site-packages/rasterio/env.py", line 401, in wrapper
    return f(*args, **kwds)
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.8/site-packages/rasterio/warp.py", line 344, in reproject
    _reproject(
  File "rasterio/_warp.pyx", line 517, in rasterio._warp._reproject
  File "rasterio/_warp.pyx", line 489, in rasterio._warp._reproject
  File "rasterio/_err.pyx", line 221, in rasterio._err.exc_wrap_pointer
rasterio._err.CPLE_NotSupportedError: Cannot find coordinate operations from 'ESRI:54009' to 'EPSG:4326'
"""

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/data/davidef/git_world/pypsa-earth/scripts/build_renewable_profiles.py", line 699, in <module>
    availability = cutout.availabilitymatrix(regions, excluder, **kwargs)
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.8/site-packages/atlite/gis.py", line 748, in compute_availabilitymatrix
    availability = list(pool.map(_process_func, shapes.index))
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.8/multiprocessing/pool.py", line 364, in map
    return self._map_async(func, iterable, mapstar, chunksize).get()
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.8/multiprocessing/pool.py", line 771, in get
    raise self._value
rasterio._err.CPLE_NotSupportedError: Cannot find coordinate operations from 'ESRI:54009' to 'EPSG:4326'

Your Environment

davide-f commented 1 year ago

This issue does not appear with the environment with gdal <= 3.6 (3.5.2).

davide-f commented 11 months ago

The issue relates to offshore wind time series (ac or dc) for a geometry that crosses the +-180 degrees.

euronion commented 11 months ago

So it only affects geometries crossing +/-180° ?

euronion commented 11 months ago

If you mention it works with gdal <= 3.6, which rasterio is installed with that gdal version? Might be rasterio related.

ollie-bell commented 3 weeks ago

Are there now any workarounds to building cutouts which cross +/- 180?

euronion commented 5 days ago

Not that I am aware of. Does the issue still persist for you? Or is the problem outdated?

davide-f commented 16 hours ago

This issue still persists. The workaround I found is avoid shapes that are multipolygons crossing the antimeridian. It would be best to split them in two with the current implementation.

A long-term solution may be to change the procedure how atlite works and avoid the image cropping completely. By exploiting a mapping function between the cutout and the geometry, it would be possible to sum only the numbers that lay wihin the region, without the need to perform the geometry operations. In PyPSA-Earth we have an implementation for computing the population and could be adapted for atlite. See: https://github.com/pypsa-meets-earth/pypsa-earth/blob/c1d4b993dc52d1f22a46efab983bcef5289afb95/scripts/build_shapes.py#L1114

That may solve the issue completely, but it requires effort