Search code examples
rrasterprojection

Is there an R function to change the projection of raster?

I am trying to plot crash data with a density interactive map however the raster doesn't add projection.

I have data frame looks like:

TA_name crashes
chch 13223
grey 2322
buller 123

Here is waht I have tried so far:

library(devtools)
install_github("mtennekes/oldtmaptools")
library(oldtmaptools)
# creates a smooth map
# bandwidth is a size of the kernel, default to 1/50th of the shortest side
crash_density <- smooth_map(crashes_TA, bandwidth = 0.5, smooth.raster = T)
tmap_mode("view")
library(viridis) # colour schemes
library(rasterVis) # raster analysis and plotting
plot3D(crash_density$raster, zfac = 0.5, col=viridis_pal(option = "B",direction = -1, alpha = 0.5))
library(spatstat)
# data to ppp class
crashes_pp <- crashes_TA %>% as.ppp()

# just filtering Christchurch City TA
chch <- TA %>%
  filter(TA_name == "chch")
# defining shape of observation window (owin)
# setting it to Christchurch
Window(crashes_pp) <- as.owin(chch)
crashes_pp
# heat map - sigma is a bandwidth, eps - final raster resolution
crashes.kd1 <- density(crashes_pp, sigma = 500, eps = 200)
# airbnb.kd1 <- density(airbnb_pp, sigma = bw.diggle, eps = 200)
image(crashes.kd1)
# multiplying
crashes.kd1 <- crashes.kd1*1000^2
# converting to raster
crashes.kd1 <- raster(crashes.kd1)
# adding crs
crs(crashes.kd1) <- "epsg:2193"

However I get an error of:

Error in CRS(SRS_string = x) : NA

Would this be due to an error in my data map

My plot runs just in wrong place due to not executing projection.

Solution

  • Here is a minimal, self-contained, reproducible example that creates the error.

    I think you should use "EPSG:2193" instead of "epsg:2193"

    Example data

    library(raster)
#Loading required package: sp
r <- raster()

    This fails

    crs(r) <- "epsg:2193"
#Error in sp::CRS(SRS_string = x) : NA

    It works (with a warning) with

    crs(r) <- "EPSG:2193"
#Warning message:
#In showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
#  Discarded datum New Zealand Geodetic Datum 2000 in Proj4 definition