2

I am trying to create density maps (rasters) from points using the smooth_map function from the tmaptools package. I would like to create one raster for each factor level of a variable (var) in my dataframe (df). All rasters should have the same extent, given by a shapefile (study_area). I want to stack them and subsequently plot them using levelplot.

However, all rasters returned by the function have slightly different extents. I could not figure out how to create the density maps with same extents. This is the code:

library(tmaptools) 
library(rgdal)
library(sf)
library(raster)

##reproducible example of dataframe
df<-setNames(data.frame(matrix(ncol = 3, nrow = 7)), c("longitude", "latitude", "var"))

df$longitude<-c(-53.30002, -50.47749, -55.85561, -57.88447, -55.83864, -58.84610, -57.49215)
df$latitude<-c(-13.037530,  -13.023480, -13.416200, -13.659120, -12.758670, -13.114460, -14.622520)
df$var<-c("A", "A", "A", "A", "B", "B", "B")

## convert var to factor
df$var<-as.factor(df$var)

##read shapefile of Mato Grosso, Brazil 
## shapefile can be downloaded from e.g. https://geo.nyu.edu/catalog/stanford-vw409fv3488
study_area<-readOGR("shape_MT.shp", layer="shape_MT") #will load the shapefile 

#create empty stack
s <- stack()

#loop over factor levels
for (i in levels(df$var)) { 
  a<-subset(df, df$var==i)
  my.sf.point <- st_as_sf(x = a, coords = c("longitude", "latitude"),crs = "+proj=longlat  +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
  density_map <- smooth_map(my.sf.point, cover = study_area) 
  r<-density_map$raster
  print(extent(r))
  s <- stack(s,r)
}

I get the following error message:

Error in compareRaster(x) : different extent

Any idea on how to solve this?

Ben
  • 28,684
  • 5
  • 23
  • 45
albren
  • 109
  • 1
  • 1
  • 8
  • Could you add a sample of data `df` to reproduce? – Ben Sep 29 '19 at 19:06
  • @Ben a reproducible example was added – albren Sep 29 '19 at 19:52
  • This is amjorly because of the difference of cell resolution between the two raster data that somehow affects the extent too (cells arenot squared). You just need to resample them before stacking. Easy fix it is if that is what you want! – Majid Oct 01 '19 at 13:50
  • Plus, as a suggestion, `+longlat` as `crs` is not the best for kernel density. – Majid Oct 01 '19 at 13:52
  • @Majid what would you suggest as crs? I tried to use UTM, but my study area includes more than one UTM zone – albren Oct 01 '19 at 15:12

1 Answers1

1

Perhaps this works for you:

## convert var to factor
df$var<-as.factor(df$var)

## read shapefile of Mato Grosso, Brazil 
## shapefile can be downloaded from e.g. https://geo.nyu.edu/catalog/stanford-vw409fv3488
study_area<-readOGR("....", layer="...")

It will also work with your own crs but perhaps ESRI:102033 South_America_Albers_Equal_Area_Conic looks like a better option here (or something similar):

crs.laea <- CRS("+proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0 +y_0=0 +ellps=aust_SA +units=m +no_defs")
study_area <- spTransform(study_area, CRS = crs.laea)

#create empty stack
s <- list()

#loop over factor levels
for (i in levels(df$var)) { 
  a<-subset(df, df$var==i)
  my.sf.point <- st_as_sf(x = a, coords = c("longitude", "latitude"), crs = "+proj=longlat +ellps=aust_SA +no_defs")

  #reprojecting to ESRI:102033 South_America_Albers_Equal_Area_Conic
  my.sf.point <- st_transform(my.sf.point, 102033)
  density_map <- smooth_map(my.sf.point, cover = study_area, nrow = 514, ncol = 510) 
  r<-density_map$raster

  extent(r) <- c(-175916.2, 1048374, 1610251, 2845008) # extent of the study area
  plot(r)

  s[[i]] <- r
}

s[["A"]] <- resample(s[["A"]],s[["B"]],method='bilinear')

stack(s)
Majid
  • 1,836
  • 9
  • 19