[R-sig-Geo] Extract Raster Data Values to Shapefile

[Julien Rodriguez]

Thanks: it is true that handling missing values when extracting from a
raster can sometimes be perilous.

Same example with missing values and playing with terra::extract arguments
to take the best decision about them.

library(terra)
library(sf)

# Create a SpatRaster
r <- rast(ncols=5, nrows=5, xmin=0, xmax=5, ymin=0, ymax=5)
val <- 1:25
# add randomly selected NAs
set.seed(220512)
val[sample(1:25,4)] <- NA

values(r) <- val

# Create a polygon object as sf
p1 <- rbind(c(1, 1), c(3, 1), c(3, 3), c(1, 3), c(1, 1))
p2 <- rbind(c(4, 0), c(4, 1), c(5, 1), c(5, 0), c(4, 0))
p3 <- rbind(c(1.6, 1.6), c(3, 1.6), c(3, 3), c(1.6, 3), c(1.6, 1.6))
myPol <- st_sf(code = c(1:3), geometry = st_sfc(list(st_polygon(list(p1)),
st_polygon(list(p2)), st_polygon(list(p3)))))

# Polygon object as SpatVector
myPolAsSpatVect <- vect(st_geometry(myPol))

plot(r)
plot(myPolAsSpatVect, add = TRUE)

# Here is your extraction, using fun argument, you can decide how to
aggregate the values by district border
terra::extract(r, myPolAsSpatVect, fun = mean)
terra::extract(r, myPolAsSpatVect, fun = mean, na.rm = TRUE)
# By default, the extraction will concern pixels whose center point is
within the polygon
terra::extract(r, myPolAsSpatVect, fun = mean, na.rm = TRUE, touches = TRUE)
terra::extract(r, myPolAsSpatVect, fun = mean, na.rm = TRUE, exact = TRUE)

# without aggregation
terra::extract(r, myPolAsSpatVect)
terra::extract(r, myPolAsSpatVect, exact = TRUE)
terra::extract(r, myPolAsSpatVect, touches = TRUE)

Regards



Le jeu. 12 mai 2022 à 09:29, Jordan Chamberlin <jordan.chamberlin using gmail.com>
a écrit :

> One way to deal with NAs is simply to add an argument to terra::extract:
> na.rm=TRUE
>
>
>
> On Wed, May 11, 2022 at 9:59 PM GilbertoCamara <
> gilberto.camara.inpe using gmail.com> wrote:
>
>> Dear Julian,
>>
>> Your script works nicely if the raster data if well-behaved. In real
>> life, raster data can have missing values, which have to be dealt with.
>>
>> Best
>> Gilberto
>> ============================
>> Prof Dr Gilberto Camara
>> Senior Researcher
>> National Institute for Space Research (INPE), Brazil
>>
>> https://urldefense.com/v3/__https://gilbertocamara.org/__;!!HXCxUKc!xb3Skg9BEorRUu-1yoW_jHkQRP9y6tC81a3GXDuYrM-JEP1U_aBOybnpLw7uroqGQsB8yGdvdj8X4xbotnvel5LZQAkw$
>> =============================
>>
>>
>>
>>
>> > On 11 May 2022, at 14:52, Julien Rodriguez <
>> julienrodriguez.ardaspm using gmail.com> wrote:
>> >
>> > Hi Luca and Zivan,
>> >
>> >
>> > Maybe this toy example using terra and sf (probably the class of your
>> > district borders object) may help.
>> >
>> > library(terra)
>> > library(sf)
>> >
>> > # Create a SpatRaster
>> > r <- rast(ncols=5, nrows=5, xmin=0, xmax=5, ymin=0, ymax=5)
>> > values(r) <- 1:25
>> >
>> > # Create a polygon object as sf
>> > p1 <- rbind(c(1, 1), c(3, 1), c(3, 3), c(1, 3), c(1, 1))
>> > p2 <- rbind(c(4, 0), c(4, 1), c(5, 1), c(5, 0), c(4, 0))
>> > myPol <- st_sf(code = c(1:2), geometry =
>> st_sfc(list(st_polygon(list(p1)),
>> > st_polygon(list(p2)))))
>> >
>> > # Polygon object as SpatVector
>> > myPolAsSpatVect <- vect(st_geometry(myPol))
>> >
>> > plot(r)
>> > plot(myPolAsSpatVect, add = TRUE)
>> >
>> > # Here is your extraction, using fun argument, you can decide how to
>> > aggregate the values by district border
>> > terra::extract(r, myPolAsSpatVect, fun = mean)
>> >
>> > Regards
>> >
>> > Julien
>> >
>> >
>> >
>> >
>> > Le mer. 11 mai 2022 à 18:19, Zivan Karaman <zivan.karaman using gmail.com> a
>> > écrit :
>> >
>> >> Hi Luca,
>> >> You should perhaps have a look at the crop and mask function in the
>> raster
>> >> or terra package.
>> >> Hope this helps!
>> >> Zivan
>> >>
>> >> On Wed, May 11, 2022 at 6:09 PM <lucfrank using mail.uni-mannheim.de> wrote:
>> >>
>> >>> Hello Everyone,
>> >>>
>> >>> I have a shapefile with district borders and a raster data file
>> (.tif).
>> >>>
>> >>> How can I extract the value from the raster data onto each district in
>> >>> the shapefile?
>> >>>
>> >>> I have the exact thing I want with two shapefiles with the command:
>> >>>
>> >>> x <- st_join(x,y,join=st_overlaps,left=TRUE,largest=TRUE)
>> >>>
>> >>> This gives me the result I want. I tried to convert the .tif file to a
>> >>> shapefile, but this didn't work for me with multiple different
>> >>> solutions I found online.
>> >>>
>> >>> I would appreciate any ideas you might have.
>> >>>
>> >>> Luca
>> >>>
>> >>> _______________________________________________
>> >>> R-sig-Geo mailing list
>> >>> R-sig-Geo using r-project.org
>> >>>
>> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!!HXCxUKc!xb3Skg9BEorRUu-1yoW_jHkQRP9y6tC81a3GXDuYrM-JEP1U_aBOybnpLw7uroqGQsB8yGdvdj8X4xbotnvel9iqWVYy$
>> >>>
>> >>
>> >>        [[alternative HTML version deleted]]
>> >>
>> >> _______________________________________________
>> >> R-sig-Geo mailing list
>> >> R-sig-Geo using r-project.org
>> >>
>> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!!HXCxUKc!xb3Skg9BEorRUu-1yoW_jHkQRP9y6tC81a3GXDuYrM-JEP1U_aBOybnpLw7uroqGQsB8yGdvdj8X4xbotnvel9iqWVYy$
>> >>
>> >
>> >       [[alternative HTML version deleted]]
>> >
>> > _______________________________________________
>> > R-sig-Geo mailing list
>> > R-sig-Geo using r-project.org
>> >
>> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!!HXCxUKc!xb3Skg9BEorRUu-1yoW_jHkQRP9y6tC81a3GXDuYrM-JEP1U_aBOybnpLw7uroqGQsB8yGdvdj8X4xbotnvel9iqWVYy$
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo using r-project.org
>>
>> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!!HXCxUKc!xb3Skg9BEorRUu-1yoW_jHkQRP9y6tC81a3GXDuYrM-JEP1U_aBOybnpLw7uroqGQsB8yGdvdj8X4xbotnvel9iqWVYy$
>>
>

	[[alternative HTML version deleted]]