0

I would like to rasterize data on protected areas for certain countries (e.g. https://www.protectedplanet.net/country/TH, download button is at the bottom of the page). The dataset contains two shapefiles, one with a multipolygon geometry and the other one with a multipoint geometry. I would like to rasterize both. While rasterizing the multipolygon geometry is straightforward, the standard functions (raster::rasterize and fasterize::fasterize) don't work with multipoint features.

How can I either rasterize the multipoint feature collection directly or convert it to something that can be rasterized with the standard functions?

2 Answers 2

2

The problem seems to be the point data's shapefile rather than the functions that seem to work with polygons and not points. You should check with the publisher of your data. I was not even able to load the point data in R. Instead, by exporting the shapefile as .csv in QGIS and loading it in R (get the csv here) as a multipoint feature, one would get:

library(rgdal)
library(raster)
# raster that will become the 'rasterized' points
r <- raster()
res(r) <- 10 # setting raster's resolution
# loading csv and making coordinates numerical:
points <- read.csv2('WDPA_Sep2019_THA-shapefile-points.csv', sep=',', stringsAsFactors = F)
points[,c('xcoord', 'ycoord')] <- data.frame(sapply(points[,c('xcoord', 'ycoord')], as.numeric))
# transforming it into a multipoint spatial feature
points <- SpatialPoints(coords = as.matrix(points[,c('xcoord', 'ycoord')]), proj4string = r@crs)

# adjusting the extent of the raster to the points'
r <- setExtent(r, points@bbox)
# rasterizing it:
r <- rasterize(points, r)

# for visualization only: set overlapping points as 1, zero otherwise:
r[!is.na(r)] <- 1
r[is.na(r)] <- 0

plot(r)
plot(points, add=T)

the outcome is the image attached. Important: I set raster's resolution to make the cells visible in the plot together with the points. That is not a needed passage.

enter image description here

2
  • Thanks! Did you try to read the shapefile with readOGR? That also didn't work in my case, but using the sf package and st_read did the job.
    – tho_mi
    Commented Sep 14, 2019 at 19:17
  • Next time try to provide the code you are working with so that we can address better your issue. Commented Sep 16, 2019 at 16:30
0

Turns out, it's not really a multipoint feature. While one can't read the file with readOGR, it works with st_read from the sf package. From there it's possible to transform the object into a SpatialMultiPoints object (shape_points below).

But looking at the features one sees that they're actually just simple points, not multipoints. The following command creates a SpatialPoints object based on the SpatialMultiPoints object:

shape_points = SpatialPoints(matrix(unlist(shape_points@coords), ncol = 2, byrow = TRUE), proj4string = CRS(projection(raster_pop)))

After that one has a SpatialPoints object which can be rasterized using the rasterize function from the raster package.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.