Can I save the weightmap class as an xarray, so I can have access to the weights using indexing?

[fllamorim]

I would like to save the weightmap class as xarray dataarray, so I can use the indexing to access the percentage of grid covered by a shapefile.
Is it possible?

Thank you in advance.

[ks905383]

Hi @fllamorim, thanks for your question.

There's no "official" way to do so (yet), but there is a workaround that I think will provide what you need.

If you run xagg with impl='dot_product' (see the docs here), then the weightmap object has an attribute overlap_da that returns a DataArray with dimensions loc (giving the pixel) and poly_idx, giving the polygon. The values of the dataarray are the relative weight given to each pixel for a given polygon. Here's what that looks like with the default data from the docs:

To add the other metadata into that, you could do the following:

overlaps = weightmap.overlap_da
for v in [col for col in weightmap.agg.columns if col not in ['pix_idxs','rel_area','coords','poly_idx']]:
    overlaps[v] = xr.DataArray(weightmap.agg[v],dims=['poly_idx'],coords = {'poly_idx':weightmap.agg['poly_idx']})

If there's interest in making this a more 'official' feature, I recommend opening an issue and encourage PRs :)

[fllamorim]

Dear Kevin,

Thank you very much for your answer. It helped me to advance a little.
I am still having problems to understand the coverage percentage values between a grid cell and the shapefile.
Sorry, it is a stupid question, but the relative ares is relative to what?

This is my grid and shapefile.

And these are the overlap values:

Why do I have no ones (1) in my coverage percentage?

If you do not mind, I am sharing my grid and shapefile with you.

grid_shapefile.zip

Sorry, another question is: do you have any idea how to reconstruct the dependent lat lon matrix from the rel_area 'loc'? This would be exactly what I need.

xarray.DataArray
'rel_area' (poly_idx: 1 loc: 238)

Coordinates:

• poly_idx (poly_idx) int64 0
• loc (loc) object MultiIndex
• lat (loc) float64 -41.25 -41.25 -41.25 -41.25 ... -33.25 -33.25 -33.25
• lon (loc) float64 -63.25 -62.75 -62.25 -61.75 ... -57.75 -57.25 -56.75

I am really sorry if this is a lot, but your library was the closest I found to solve my problem.

Best regards,

Felipe

[ks905383]

Hi Felipe,

Sorry, it is a stupid question, but the relative ares is relative to what?

Relative areas are calculated from the perspective of the polygon. So for each pixel up there, the relative area shows how much the overlap between that pixel and the polygon makes up of the whole polygon's area. (If you summed up the weights up there, they would add up to 1)

To visualize that a bit more clearly, there's a way to create a diagnostic figure that shows this explicitly. Here's what that would look like for your data:

(the diagnostic figure normalizes even further by the relative area / weight of the pixel with the largest overlap, but this is essentially what that overlap_da is showing.)

Sorry, another question is: do you have any idea how to reconstruct the dependent lat lon matrix from the rel_area 'loc'? This would be exactly what I need.

You could always do .unstack() but that would I think just get you the same as the overlap_da above.

I would like to save the weightmap class as xarray dataarray, so I can use the indexing to access the percentage of grid covered by a shapefile.

Ah, I think I now understand what you mean by this. You want to know what % of the grid is covered by the polygon? (May I ask what you would use that number for?)

In that case, you could use xa.create_raster_polygons() to create polygons for each pixel, and then use gpd.overlay() from geopandas to get the intersection of the two geodataframes.

# Get polygons from every grid cell
raster_polys = xa.core.create_raster_polygons(ds)

# Get intersection between polygon and the grids, creating a geodataframe
# where each row is a polygon representing the intersection between a 
# grid cell and the polygon
poly_intersect = gpd.overlay(gdf,raster_polys['gdf_pixels'])

# Get percentage of grid covered by the polygon, by 
# - changing to a crs with units of km (the default crs used here has units of degrees)
# - calculating the area of each row of the intersection / the raster grid
# - summing over both and getting their ratio
tot_rel_area = (poly_intersect.to_crs('EPSG:3857').area.sum() / 
                           raster_polys['gdf_pixels'].to_crs('EPSG:3857').area.sum())

[fllamorim]

Dear Kevin,

Thank you very much. With your help I could solve my problem.

May I ask what you would use that number for?

This is a sensitive analysis I am doing for absence and pseudo absence related to chosen random grid cells and polygons. I need the covered percentage of a grid cell by a polygon just to normalize the analysis.

Thank you very much again for all your time.

Best regards,

Felipe