Add handling of AOI as computational region, GRASS map or GeoJSON

[ninsbl]

from pathlib import Path

from shapely import from_wkb
from shapely.geometry import shape
# from shapely.geometry.polygon import Polygon

import grass.script as gs

from grass.pygrass.vector import VectorTopo

def get_aoi(aoi):
    """Get the AOI for querying"""
    # Handle empty AOI
    if not aoi:
        reg = gs.parse_command("g.region", flags="gl")
        return tuple(map(float, [reg['sw_long'], reg['sw_lat'], reg['ne_long'], reg['ne_lat']]))

    # Handle AOI as GeoJSON (GeoJSON is EPSG:4326 by definition)
    # A solution with GDAL/OGR would be a (probably more flexible) alternative
    # However GeoJSON is probably more straight forward to use as it is WGS84 by standard definition
    aoi_path = Path(aoi)
    if aoi_path.exists():
        try:
            geo_json = json.loads(aoi_path.read_text())
            # Todo: handle multiple geometries (use firt give warning
            #       maybe handle geometry types (Multi, line, point, ...)
            return shape(geo_json)
        except RuntimeError:
            gs.fatal(_("Cannot read input file <{}>").format(aoi))

    # Handle AOI as GRASS vector map
    else:
        aoi_map = VectorTopo(aoi)
        try:
            aoi_map.open("r")
            aoi_geoms = vmap.areas_to_wkb_list()
            if len(aoi_geoms) > 1:
                gs.warning(_(
                    "Found more than one AOI geomtry in vector map <{}>. Using only the first.
                ).format(aoi))
        except RuntimeError:
            gs.fatal(_("Could nor read AOI from vector map <{}>).format(aoi))
        aoi_geom = aoi_geoms[0][-1]
        # Here we need to reproject if location CRS != EPSG:4326
        # Note that shapely geometries are unaware of the CRS
        return from_wkb(aoi_geom)

Note that the code is not tested...

[veroandreo]

What about using GRASS raster or vector maps as aoi too?

[HamedElgizery]

@ninsbl

# Here we need to reproject if location CRS != EPSG:4326
# Note that shapely geometries are unaware of the CRS

I think the way to do that would be to:

1. Create a temp_location with EPSG:4326 CRS.
2. Reproject the targeted vector into the temp_location.
3. Open the new reprojected vector, and get the wkb value from it.
4. Delete the temp_location.
5. Create the shapely file with the WKB and return it.

I don't think that is really clean, but I couldn't find a way to reproject the Vector in Python on the go... Can you refere me to any references, if any, that could help reach a cleaner approach, or if this is good enough I will implement it...

[ninsbl]

I would use something along these lines:
https://pcjericks.github.io/py-gdalogr-cookbook/projection.html

It adds GDAL (osr/ogr) as an external dependency, but that should be OK.
You can get the source crs with g.proj.

[HamedElgizery]

Thanks. I have been trying to use it, but there is some trouble because Shapely uses (longitude, latitude), while GDAL uses (latitude, longitude)... We could easily convert one to the other, but the problem is that when I try to reproject a vector from EPSG:4326 to EPSG:3358, I do have to convert the shapely object from (longitude, latitude) to (latitude, longitude), but when I try to reproject the other way around (from EPSG:3358 to EPSG:4326), I don't need to.... So apart from this case, I think that there might be other edge cases, that might cause other problems... (Also there is some slight errors with the coordinates, but I think it just a precision problem).

Also, I am still not sure why I don't need to convert the coordinates order, when converting from EPSG:3358 to EPSG:4326.

What about we use a similar implementation like in i.sentinel.download.
https://github.com/OSGeo/grass-addons/blob/grass8/src/imagery/i.sentinel/i.sentinel.download/i.sentinel.download.py#L246

Or if this doesn't have a major drawbacks?

• Create a temp_location with EPSG:4326 CRS.
• Reproject the targeted vector into the temp_location.
• Open the new reprojected vector, and get the wkb value from it.
• Delete the temp_location.
• Create the shapely file with the WKB and return it.

[veroandreo]

IMO, something on the lines i.sentinel.download does could work, in that way we do not need any extra dependency nor creating a temp-location. What do you think, @ninsbl and @lucadelu?

[HamedElgizery]

Implemented it in a very similar way to i.sentinel.download method. OSGeo#1097

[HamedElgizery]

@ninsbl Regarding the GeoJSON file support, during the last meeting we discussed that if a user have a GeoJSON file, it can first be imported as a vector map, and that vector map could then be used to get the AOI, and so that would be a way to use GeoJSON indirectly.
Do you think that it is enough to have it that way, or would it be better if direct GeoJSON support was added as well?

[HamedElgizery]

Closing this issue.