Skip to content

Polygon Data into Grid Structure

3 messages · iucir@@k m@iii@g oii m@ii@u@i-m@@@heim@de, Roger Bivand, Michael Sumner

Original
iucir@@k m@iii@g oii m@ii@u@i-m@@@heim@de
# 
Hello everyone,

I have the Aqueduct 3.0 Data from the WRI Institut which is in  
Multipolygon format.

To perform a regression analysis with different datasets I want to  
transform these polygons into a grid.
Basically, I want to overlay a grid over the map and assign each cell  
the value of the underlying polygon.
In a second step, I want to count how many points (different dataset)  
are inside each cell.


I have tried with the sp commands GridTopology, and SpatialPoints. I  
can produce the grid but fail to crop with my map.
Also, it is not clear to me how I can assign the underlying value of  
the polygon.

This is the code I have tried so far:
(I have copied this from doing a Course on DataCamp which I can't  
completely copy exactly because my data is different)

# -25,-49  my lower left corner (these are coordinates)
grid <- GridTopology(c(-25,-49),c(1, 1), c(72, 48))


gridpoints <- SpatialPoints(grid, proj4string = CRS(projection(aqua)))
plot(gridpoints)

cropped_gridpoints <- crop(gridpoints, aqua)
plot(cropped_gridpoints)

spgrid <- SpatialPixels(cropped_gridpoints)
plot(spgrid)


I hope someone has an idea to help me solve this.
Thanks to everyone who has any insights into this

Luca Frank
Roger Bivand
#
On Mon, 2 May 2022, lucfrank at mail.uni-mannheim.de wrote:

            

      
      
      
    

        
This is not clear, neither aqueduct nor WRI mean anything to me and 
possibly most others.

            

      
      
      
    

        
So a standard intersection between a raster grid and polygons - but you 
are modifying the support of the data quite brutally. If you need to 
integrate data with varying support, do not treat the estimated values for 
the new entities as fixed, they are estimates with error distributions.

Why not simply count the points within each polygon?

            

      
      
      
    

        
Not a robust idea; build up your analysis from the support of your own 
data based on robust sources.

            

      
      
      
    

        
This looks rather odd, seeming to use geographical coordinates, so the 
areas in the grid cells vary with latitude (as would the counts).

            

      
      
      
    

        
crop is from terra (in older workflows raster), not sp. Perhaps using 
terra for the whole exercise is advisable? Your workflow template doesn't 
seem to be current, or even adequate, I'm afraid.

Roger

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Michael Sumner
    

    

      
      #
    

  


  

      
On Tue, May 3, 2022 at 2:22 AM <lucfrank at mail.uni-mannheim.de> wrote:

        

            

      
      
      
    

        
try

library(terra)
p <- vect(<yourpolyfile>)
r <- rasterize(p, rast(<your-raster-file>), field = seq_along(p))

In a second step, I want to count how many points (different dataset)

            

      
      
      
    

        
I'd use terra, and terra::cells(grid, vect(points))

where 'grid' is from

rast(<your-raster-file>)

and points is cbind(x, y) just your points. If grid and points are in
different projections you'll need to declare that for the points in vect( ,
crs  = ) unless you read it from a source with that info.

You get the cell id of the raster, which you can store alongside the x-y
values and summarize from there.

Best, Mike