0

I have access to a tif file which contains satellite images of the Netherlands (e.g. https://ns_hwh.fundaments.nl/hwh-ortho/2022/Ortho/1/01/beelden_tif_tegels/2022_079000_449000_RGB_hrl.tif), and a shapefile containing the plot area (parcel area) of houses.

area <- readOGR( 
  dsn = paste0( getwd(), "/Perseel/" ) , 
  layer = "Essellanden_percelen",
  verbose = F )

plot(area)

The .shp file looks like the following: Plot area

My end goal is to output each polygon of the plot area separately as an image (.jpg or .png) to be able to use them in image recognition, but I don't know how to get there.

The first step would be to overlay the plot areas on top of the satellite image, such as: Example map

Next step would be to cut out the plot areas one by one and output them as .jpg.

I read the data into R as a RasterBrick (.tif) and SpatialPolygonsDataFrame (.shp) and tried to convert either to be of the same data type so it would be easier to lay on top of each other, but I was unable to do this.

M--
  • 25,431
  • 8
  • 61
  • 93
Maaike
  • 11
  • 2

1 Answers1

0

Assuming that your tif raster is georeferenced, you can use raster::crop to cut the image to the extent of each of the polygons in your shapefile by looping over it.

library(rgdal)
library(raster)

lapply(seq_len(nrow(area)) function(i) crop(aerial_image, extent(area[i,), snap="out")
M--
  • 25,431
  • 8
  • 61
  • 93
  • Thank you for the suggestion! I tried your code and adjusted it as follows: ``` loop <- function(i) { crop(image, extent(area[i,], snap = "out" ) ) } lapply( seq_len ( nrow( area ) ), loop) ``` This is giving me an error: "Cannot get an Extent object from argument y". Running only "crop(image, extent(area[1,]))" gives the error that the extents do not overlap. Do you know how I can check the extent and make sure they are similar? – Maaike Apr 11 '23 at 15:02
  • In addition the projection of area is: "+proj=longlat +datum=WGS84 +no_defs". The projection of aerial_image is: "+proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs" – Maaike Apr 11 '23 at 15:09
  • @Maaike set the projections to be the same: https://stackoverflow.com/questions/50232331/set-the-right-crs-on-sf-object-to-plot-coordinate-points Moreover, no need to save the clip as a function. Use it as is in my answer, it will loop through your data. I am not sure what did you do with `loop` and why did you do it. After you set the projections, it should work. Based on your screenshot, the extents are matching. – M-- Apr 11 '23 at 16:00