Intersection of polygons with raster

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Intersection of polygons with raster

Kent Johnson
Hi,

I have a raster object `flood` containing projected flooding levels and a
simple features object `parcels` containing MULTIPOLYGONs representing
property parcels. I would like to find all the parcel polygons for which
there is any flood > 1 foot. What function(s) can do this? Something like
raster::intersect(parcels, flood >= 1)

I'll put together a reprex if that helps. Hoping someone can point me to
the right function to do this, I haven't found anything promising.

Thanks,
Kent

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Intersection of polygons with raster

Michael Sumner-2
Hi Kent, this is pretty straightforward with raster::extract. If speed is
an issue you can burn the polygon id into the grid and then group values by
that.

library(raster)
r <- raster(volcano)

## very simplistic polygon layer example
r2 <- raster(r)
res(r2) <- res(r2) * 20
p <- sf::st_as_sf(as(r2, "SpatialPolygonsDataFrame"))

## a dummy threshold value
onefoot <- 150
## grab the first column from extract (it gives multiple columns for
multi-layer rasters)
p$onefoot <- extract(r, p, fun = function(x, na.rm = TRUE) any(x >
onefoot))[,1]

plot(sf::st_geometry(p), col = p$onefoot + 1)
contour(r, levels = onefoot, col = "white", add = T)


## if speed is an issue, use fasterize for a more abstract workflow
library(fasterize)
p$rownum <- 1:nrow(p)
idgrid <- fasterize(p, r, field = "rownum")
p$onefoot <- tapply(values(r), values(idgrid), function(x, na.rm = TRUE)
any(x > onefoot))
plot(p)

FWIW it's always useful to provide a reprex, since that does take time and
might miss the mark for your situation. It's also trickier to provide two
objects that are relevant to each other, so any upfront work you can do in
an example helps the answerer.

HTH

On Sat, 19 May 2018 at 12:20 Kent Johnson <[hidden email]> wrote:

> Hi,
>
> I have a raster object `flood` containing projected flooding levels and a
> simple features object `parcels` containing MULTIPOLYGONs representing
> property parcels. I would like to find all the parcel polygons for which
> there is any flood > 1 foot. What function(s) can do this? Something like
> raster::intersect(parcels, flood >= 1)
>
> I'll put together a reprex if that helps. Hoping someone can point me to
> the right function to do this, I haven't found anything promising.
>
> Thanks,
> Kent
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Intersection of polygons with raster

Kent Johnson
Thank you Michael, that is exactly what I needed. It is fairly slow for my
case but I only need to do it once :-)

And thank you for providing an example, I'm not too familiar with rasters
and would have had some trouble creating it.

Kent

On Fri, May 18, 2018 at 11:41 PM, Michael Sumner <[hidden email]> wrote:

> Hi Kent, this is pretty straightforward with raster::extract. If speed is
> an issue you can burn the polygon id into the grid and then group values by
> that.
>
> library(raster)
> r <- raster(volcano)
>
> ## very simplistic polygon layer example
> r2 <- raster(r)
> res(r2) <- res(r2) * 20
> p <- sf::st_as_sf(as(r2, "SpatialPolygonsDataFrame"))
>
> ## a dummy threshold value
> onefoot <- 150
> ## grab the first column from extract (it gives multiple columns for
> multi-layer rasters)
> p$onefoot <- extract(r, p, fun = function(x, na.rm = TRUE) any(x >
> onefoot))[,1]
>
> plot(sf::st_geometry(p), col = p$onefoot + 1)
> contour(r, levels = onefoot, col = "white", add = T)
>
>
> ## if speed is an issue, use fasterize for a more abstract workflow
> library(fasterize)
> p$rownum <- 1:nrow(p)
> idgrid <- fasterize(p, r, field = "rownum")
> p$onefoot <- tapply(values(r), values(idgrid), function(x, na.rm = TRUE)
> any(x > onefoot))
> plot(p)
>
> FWIW it's always useful to provide a reprex, since that does take time and
> might miss the mark for your situation. It's also trickier to provide two
> objects that are relevant to each other, so any upfront work you can do in
> an example helps the answerer.
>
> HTH
>
> On Sat, 19 May 2018 at 12:20 Kent Johnson <[hidden email]> wrote:
>
>> Hi,
>>
>> I have a raster object `flood` containing projected flooding levels and a
>> simple features object `parcels` containing MULTIPOLYGONs representing
>> property parcels. I would like to find all the parcel polygons for which
>> there is any flood > 1 foot. What function(s) can do this? Something like
>> raster::intersect(parcels, flood >= 1)
>>
>> I'll put together a reprex if that helps. Hoping someone can point me to
>> the right function to do this, I haven't found anything promising.
>>
>> Thanks,
>> Kent
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> <https://maps.google.com/?q=203+Channel+Highway+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
> Kingston Tasmania 7050 Australia
> <https://maps.google.com/?q=203+Channel+Highway+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo