Search code examples
rr-rasterr-leafletr-mapview

Raster output from rasterized sf points does not align with points

I am trying to rasterize some points and am getting a mismatch between the points and the rasters despite the crs being the same. If I convert the raster to polygons it lines up perfectly with the sf points data, but I can't figure out why the raster doesn't.

library(spData)
library(sf)
library(raster)
library(mapview)

## import some data
cycle_hire_osm = spData::cycle_hire_osm

## project to metres
cycle_hire_osm_projected = st_transform(cycle_hire_osm, crs = 27700)

## create raster template to rasterize to
raster_template <- raster(extent(cycle_hire_osm_projected), nrows = 10, ncols = 10, crs = 27700)

## rasterize the points
ch_raster1 = rasterize(cycle_hire_osm_projected, raster_template, field = 'capacity',
                       fun = sum, crs = 27700)

## convert raster to polygons
ch_poly <- rasterToPolygons(ch_raster1)

If these are plotted there are raster cells that have a value but have no points in.

## plot on a map
mapview(ch_poly)+cycle_hire_osm_projected+ch_raster1

Additional example based on reply to show the output as base, mapview and leaflet (note: I had to install the development versions of mapview and leaflet in order to plot SpatRasts)

library(spData)
library(sf)
library(terra)
library(dplyr)

# remove NAs so they are not considered
dat <- spData::cycle_hire_osm %>% filter(!is.na(capacity))
v <- vect(dat)

r <- rast(v, nrows=10, ncols=10)

chr <- rasterize(v, r, field="capacity", fun=sum, na.rm=TRUE)

## base plot
plot(chr)
points(v, cex=.5)
points(v[is.na(v$capacity)], cex=.5, col="red")



## mapview
library(mapview)
mapview(v)+chr

enter image description here

##leaflet
library(leaflet)

leaflet() |> 
  addProviderTiles(providers$CartoDB.Positron) |>
  addCircles(data = v) |>
  addRasterImage(chr)

All three plots are the same raster but the raster appears to have a different number of cells with values in each plot?

Adding example with project = FALSE as explain by @RobertHijmans

leaflet() |> 
  addProviderTiles(providers$CartoDB.Positron) |>
  addCircles(data = v) |>
  addRasterImage(chr, project = FALSE)

enter image description here

Solution

  • I do not see that issue with "raster" nor with its replacement, "terra" when using base-plot

    library(spData)
library(terra)
# using a SpatVector for easier plotting; results are the same
v <- vect(spData::cycle_hire_osm)
v <- project(v, "epsg:27700")
r <- rast(v, nrows=10, ncols=10)

chr <- rasterize(v, r, field="capacity", fun=sum, na.rm=TRUE)

plot(chr)
points(v, cex=.5)
points(v[is.na(v$capacity)], cex=.5, col="red")

    enter image description here

    But note that there are some cells with values of zero where all values of v$capacity are NA. That is because

    sum(NA, na.rm=TRUE)
#[1] 0

    To avoid that from happening you could do

    vv <- v[!is.na(v$capacity)]
chr <- rasterize(vv, r, field="capacity", fun=sum, na.rm=TRUE)

    The reason you see differences when using mapview/leaflet is that these use transform your data to the crs that they use. To avoid that use the Pseudo-Mercator (EPSG:3857) crs, and, in leaflet, use project=FALSE when adding the raster data.

    addRasterImage(chr, project=FALSE)