Is there a way to get a list of cells that intersect with a polygon-agent?

[ifoxfoot]

Hello!

My goal is to add elevation to a raster layer in every cell that intersects with an agent. The agent is a large polygon so it intersects with lots of cells.

Here is an non-functioning code example within the step function of my polygon-agent

        if self.model.step_count%4 == 0:
            #get shell weight 
            shell_weights = [a.shell_weight for a in self.model.space.get_intersecting_agents(self) 
                                       if isinstance(a, (Oyster, Shell))]
            self.total_shell_weight = sum(shell_weights)
            #amount to raise reef by
            self.mm_of_growth = (((self.total_shell_weight * 100000)/1000)/self.SHAPE_Area)/0.731
            #Get the cells that intersect with the polygon-agent
            self.model.space.raster_layer.get_intersecting_cells(self).new_elevation =+ self.mm_of_growth

I'm stuck on the last line, but basically I just need to get every cell that intersects with the agent. The get_intersecting_cells(self) part was wishful thinking...but there must be a way to do this?

Thank you!

[wang-boyu]

I can think of two ways of doing it at the moment:

1. Use the mask operation from rasterio, see Masking a raster using a shapefile. But I'm not sure how you can mask a polygon directly on top of a numpy array, which is what we have in Mesa-Geo. As a workaround, you may want to:

   • open up the raster layer using rasterio when initializing the model, and save a copy aside as a dummy index layer.
   • Then at each step you can reset its attribute values to mm_of_growth, and mask the polygon agent using this dummy index layer, to get the raster cells that intersect with the polygon.
   • Thereafter this could be as easy as a matrix addition - the area masked out will have attribute of 0 (if I'm not wrong), and you can add this masked matrix all at once to your raster layer values.

2. Another way is to convert each cell to a Point, then use vector operations to see if it's within the polygon. An example is in the population example model:

   mesa-geo/examples/population/population/model.py

   Lines 72 to 78 in c41d1ae

   	for cell in self.space.population_layer:
   	popu_round = math.ceil(cell.population)
   	if popu_round > 0:
   	for _ in range(popu_round):
   	num_agents += 1
   	point = Point(self.space.population_layer.transform * cell.indices)
   	if not point.within(self.space.lake):

   Here I'm testing whether a person in a raster cell is not in a polygon (a lake). However, this might be very slow, since you'll need to iterate through each single cell at every step.

[ifoxfoot]

Hi @wang-boyu thanks for the hint. I decided to go with option 1 (kind of) and I feel like I am soooooo close. Instead of masking the raster by the polygon-agent I chose to rasterize the shape of the polygon agent so that it has the values I want inside the polygon and zeros everywhere else. Outside of mesa, my workflow is as follows:

import numpy as np
import geopandas as gpd
from rasterio.features import rasterize
import rasterio 
from shapely.geometry import Polygon
import pandas as pd
import matplotlib.pyplot as plt


# create a geopandas dataframe from the polygons and attributes
# create a list of polygons
polygons = [Polygon([(0, 0), (0, 1), (1, 1), (1, 0)]),
            Polygon([(1, 1), (1, 2), (2, 2), (2, 1)]),
            Polygon([(2, 2), (2, 3), (3, 3), (3, 2)])]

# create a pandas dataframe with some attributes
df = pd.DataFrame({'id': [1, 2, 3],
                   'name': ['polygon 1', 'polygon 2', 'polygon 3'],
                   'value': [10, 20, 30]})

# create a geopandas dataframe from the polygons and attributes
gdf = gpd.GeoDataFrame(df, geometry=polygons)

# specify the dimensions and resolution of your output raster
width = 100
height = 100
transform = rasterio.transform.from_bounds(*gdf.total_bounds, width=width, height=height)

# create an empty raster to store the output
image = np.zeros((height, width), dtype=np.uint8)

# rasterize the polygons in the geodataframe and store the result in the output raster
shapes = ((geom, 255) for geom in gdf.geometry)
ras = rasterize(shapes=shapes, out=image, transform=transform)

# plot rasterized
plt.imshow(ras)
plt.show()

within my model, the code looks like this:

https://github.com/ifoxfoot/oyster_model/blob/5590a4fbcb7a11888b5afed0683569cf8e527c1c/agents.py#L261-L278

but when I run it I get the error

  File "/Users/irisfoxfoot/Desktop/oyster_model/agents.py", line 269, in step
    vals = ((geom, self.mm_of_growth) for geom in self.geometry)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: 'MultiPolygon' object is not iterable

which makes sense because the the reef agent is created from a shapefile with multiple polygons, and each polygon becomes a reef agent. Is there a way for the agent to select the geometry of only the polygon it is in?

[wang-boyu]

Would it work if you use self.geometry directly? That is, remove these lines:

https://github.com/ifoxfoot/oyster_model/blob/5590a4fbcb7a11888b5afed0683569cf8e527c1c/agents.py#L268-L269

And change these lines:

https://github.com/ifoxfoot/oyster_model/blob/5590a4fbcb7a11888b5afed0683569cf8e527c1c/agents.py#L274-L275

into something like

#rasterize polygon
ras = rasterize(
    shapes=self.geometry,
    fill=0,                          #  <- fill empty area outside polygon to 0
    out=self.model.space.raster_layer,
    transform=transform,
    default_value=self.mm_of_growth  #  <- set values inside polygon
)

I've never used this rasterize function before, hopefully this change does what it's supposed to.

Plus, I'm not sure whether the following lines will work, and I might be wrong:

https://github.com/ifoxfoot/oyster_model/blob/5590a4fbcb7a11888b5afed0683569cf8e527c1c/agents.py#L276-L277

You probably would need to call get_raster() from the raster layer to get the numpy array, add the new matrix, and call apply_raster() in the raster layer to update the new values.

[ifoxfoot]

@wang-boyu , I finally got it working!

I ended up having the agent initialize with its own polygon, since there are multiple agents of the same type that were each associated with a different polygon.

self.shape = poly.loc[poly['OBJECTID'] == self.unique_id]

Then, in the step function, I rasterized the polygon with the mm_of_growth value like so.

https://github.com/ifoxfoot/oyster_model/blob/f7aa37de2facfa48f0b9647dcc4604400cd5587c/agents.py#L274-L282

            vals = ((geom, self.mm_of_growth) for geom in self.shape.geometry)
            #define transform
            transform = rio.transform.from_bounds(*self.model.space.raster_layer.total_bounds, 
                                                  width=self.model.space.raster_layer.width, 
                                                  height=self.model.space.raster_layer.height)
            # create an empty raster to store the output
            image = np.zeros((self.model.space.raster_layer.height, self.model.space.raster_layer.width), dtype=np.float32) 
            #rasterize polygon
            ras = rasterize(shapes=vals, out=image, transform=transform)

Then, you were right about using get_raster and apply_raster.

https://github.com/ifoxfoot/oyster_model/blob/master/agents.py#L284-L287:~:text=ras%20%2B%20elev,%22elevation%22)

            elev = self.model.space.raster_layer.get_raster("elevation")
            new_elev = ras + elev
            self.model.space.raster_layer.apply_raster(new_elev, 
                                                       attr_name = "elevation")

and it worked! Thank you for your help :)

[wang-boyu]

Excellent!

I got a question regarding your reef agent initialization though -

I ended up having the agent initialize with its own polygon, since there are multiple agents of the same type that were each associated with a different polygon.

self.shape = poly.loc[poly['OBJECTID'] == self.unique_id]

Is this self.shape the same as the self.geometry in reef agents? If yes then I'm guessing you can access it directly through self.geometry after creating these geoagents from geodataframe.