Raster Map and on land check

[tommasogallingani]

Hi all,
I've one question about the reason why during the modeling step, the spill start position is checked on map using the on_land_pixel method in the RasterMap class (that is the father class for MapFromBNA). It seems that in that way the "on land" position could be influenced by the raster size that can be passed to the class during init.
Why the "spillability" is not checked on Polygons instead?
Thanks

    def _on_land_pixel(self, coord):
        """
        returns 1 if the point is on land, 0 otherwise

        :param coord: pixel coordinates of point of interest
        :type coord: tuple: (row, col)

        .. note:: Only used internally or for testing -- no need for external
                  API to use pixel coordinates.
        """
        # if pixel coords are negative, then off the raster,
        # so can't be on land
        if self._off_raster(coord):
            return False
        else:
            return self.raster[coord[0], coord[1]] & self.land_flag

    def on_land(self, coord):
        """
        :param coord: (long, lat, depth) location -- depth is ignored here.
        :type coord: 3-tuple of floats -- (long, lat, depth)

        :return:
         - 1 if point on land
         - 0 if not on land

        .. note:: to_pixel() converts to array of points...
        """
        return self._on_land_pixel(self.projection.to_pixel(coord,
                                                            asint=True)[0]) ```

[ChrisBarker-NOAA]

The rasterized map is used during the model run for beaching -- so we want to use exactly the same map for checking the "spillability" as for beaching.

Yes, this is influenced by the resolution of the raster -- but so is the entire beaching process, so it makes sense.

We've wanted for years to use a beaching algorithm that uses the full vector map -- but that is computationally substantially more expensive. A simple implementation would be O(N) in the number of polygon segments, and even with a smart spatial index, it would be O(logN) -- for complex polygons, that would be a substantial performance hit.

With modern computers, you can increase the resolution of the raster substantially if need be.

[tommasogallingani]

Thanks @ChrisBarker-NOAA. I was implementing a shape file based map instead of the default BNA and I fall into this issue. But now it makes sense to me.

Is it the same also for mover? Let’s say that I have a grid current mover, the evaluation of “delta” done for each step is done using raster or passed mover grid? How extrapolation is handled in that case(e.g. if the grid is not extended up to the “rasterized” coastline)?

thanks
Tommaso

[ChrisBarker-NOAA]

Ahh -- we've been meaning to add a way to build a map from a shapefile (or geojson, or ...) If you share you code, be way be able to include it -- should be pretty straightforward.

Is it the same also for mover? Let’s say that I have a grid current mover, the evaluation of “delta” done for each step is done using raster or passed mover grid?

The rasterized map is only used for beaching. the Movers do whatever they will do, unaware of the map. For a Grid Mover, the native grid is interpolated to get the velocity at the time and place of the element. If the element is off the grid, it will get a velocity of zero.

If there is a gap between the velocity grid and the land map, the element could get stuck, but usually there's a random diffusion and/or wind that will move it along.

[tommasogallingani]

Thanks for the details about mover working principle!

Ahh -- we've been meaning to add a way to build a map from a shapefile (or geojson, or ...) If you share you code, be way be able to include it -- should be pretty straightforward.

Yes of course. Now the basic implementation is to load a land polygon shape file and pass each polygons to the Polygon class as has been implemented for the BNA polygons. But we are working in order to speed up the process. I'll let you know.

Of course the Outputter Renderer will have the same implementation.

[ChrisBarker-NOAA]

Of course the Outputter Renderer will have the same implementation.

Actually, we've been meaning to refactor that -- at this point, the Map stores the land polygons anyway (set up so they can be passed to the Web client), so the Renderer should be refactored to use the polygons that the Map stores anyway, so no need to duplicate that code.

e.g. the Renderer can use the Map.get_polygons() method:

def get_polygons(self):
        polys = {}
        polys['spillable_area'] = self.spillable_area
        polys['map_bounds'] = self.map_bounds
        polys['land_polys'] = self.land_polys
        return polys

[tommasogallingani]

hi @ChrisBarker-NOAA, we are now testing a maps from polygons function that works pretty well with multi polygons shape file map like the one available here for the European Coast Line.

At the moment we create a child class (MapFromShape) that inherit method from the standard MapFromBNA class. The same should be done for Renderer.
The only difference with respect to ReadBNA functions is that the return is a tuple of Outputs Gnome Polygons (as the standard function) and the list of polygons as geopandas series. The latter can be used to eval spill ability even for polygons using the difference methods of geo polygons.

import warnings
import numpy as np
import geopandas as gpd
from gnome.maps import RasterMap, MapFromBNASchema
from gnome.maps.map import ShiftLon360, ShiftLon180, MapFromBNA
from gnome.outputters import Renderer, Outputter
from gnome.utilities import projections
from gnome.utilities.geometry.polygons import PolygonSet
from gnome.utilities.map_canvas import MapCanvas


def ReadShape(filename, polytype="list", dtype=np.float64):
    """
    Read a shape file and load as Gnome would expect polugons

    :param filename:
    :param polytype:
    :param dtype:
    :return:

    """


    if polytype == 'list':
        raise NotImplementedError()
    elif polytype == 'PolygonSet':
        from gnome.utilities.geometry import polygons
        Output = polygons.PolygonSet(dtype=dtype)


        df = gpd.read_file(
            filename).explode()
        vals = df['geometry'].values
        for idx_, val in enumerate(vals):
            poly = (np.array(list(zip(*val.exterior.coords.xy))), 'polygon', str(idx_), '1')
            Output.append(poly[0], poly[1:])

    else:
        raise ValueError('polytype must be either "list" '
                         'or "PolygonSet"')

    return Output, vals

[ChrisBarker-NOAA]

Thanks! I"ll try to give this a test soon.

However, would it be possible to refactor it to use shapely, rather than geopandas? WE already have a shapely dependency, and it would be nice to have one less.

@jay-hennen: This might be an opportunity to refactor the Map code to jsu tuse shapely for everything. When I wrote that code years ago, Shapely was a troublesome dependency (or maybe even didn't exist yet!) -- maybe it's time to just buy into it.

@tommasogallingani: As I suggested earlier, I think we should refactor the Renderer code to use the polygons from the Map object, instead of re-loading it from the vector file -- that would be one less place that has to understand file formats :-)

[tommasogallingani]

Hi Chris, I'll give it a try also with shapely. Concerning the renderer, you cannot just pass land polygons to the renderer getting the from map?

  if hasattr(getattr(model, 'map', None), 'land_polys'):
      map_file = None
      land_polygons = model.map.land_polys
  else:
      map_file = map_filepath
      land_polygons = None
  render = Renderer(
      output_timestep=gs.minutes(model.time_step / 60),
      image_size=(640, 480),
      viewport=bounds,
      formats=['png'],
      land_polygons=land_polygons,
      map_filename=map_file,
      output_dir=render_dir
  )

[jay-hennen]

I apologize for the insufficient documentation for the Renderer.

The land_polygons keyword argument is expecting an iterable (list/set) of gnome.utilities.polygons.Polygon objects. I believe the gnome implementation of the Polygon does not duck type to the shapely or geopandas Polygon (since it either predates or is contemporary with them).

In my opinion it would be an improvement if we replaced the gnome.utilities.polygons objects generally with shapely objects. The trouble of course, is it's been a refactor with very little or vague extrinsic benefit

[ChrisBarker-NOAA]

you cannot just pass land polygons to the renderer getting the from map?

Yes, you can! Sorry about that, I had forgotten we'd added that.

@jay-hennen: I think our examples and demo scripts don't take advantage of that -- we should go in and edit those!

As for Shapely polygons internally, yes no duck-typing replacement, though we may be able to add a little wrapper or something that would extract the coords from Shapely Polygons. But as you say, this is an internal detail -- I dont think it's exposed anywhere to scripting (Or save files?) so we could refactor that without breaking anything.

[ChrisBarker-NOAA]

@tommasogallingani:

I've just updated (in the develop branch) the Renderer to be able to take a GnomeMap object, rather than a BNA filename for the land polygons.

This call GnomeMap.get_polygons(), which in turn uses self.land_polys, which as Jay said, needs to be an iterable of gnome.utilities.geometry.polygons.Polygon objects (or a gnome.utilities.geometry.polygons.PolygonSet.

These are not hard to create from numpy arrays of points, so you should be able to read the shapefile, and then make these polygons, and have the whole thing integrate with the rest of the system.

Note that refactoring to use Shapely everywhere would be good, but it's a bigger refactor than we're going to do anytime soon.

Final Note: I think a "polygons_from_file" that can read BNA, shapefile (and geoJSON, and ...) could be handy -- then we could just call that everywhere we need the polygons, and turn MapFromBNA into simply MapFromVectorFile.

PRs welcome :-)