Search code examples
rleafletraster

How to set projection of interpolated raster for plotting in leaflet plot for r?


I'm trying to add a raster image to a leaflet. While not totally happy with my interpolation, I do have one that plots through the the normal plot function. The code I have shown below should be able to run on its own, as I provided sample data for the input data.

ozone_df<-data.frame("Longitude"=runif(200, -112.245075*10000,-111.455581*10000)/10000)
ozone_df$Latitude<-runif(200, 40.063614*10000,40.827281*10000)/10000
ozone_df$Ozone<-runif(200, 0,115)

#create grid tick marks
small_grid_x = seq(-111.455581,-112.245075,length.out=500)
small_grid_y = seq(40.063614,40.827281,length.out=500)
#create grid nodes
krig_grid_small<-expand.grid(small_grid_x,small_grid_y)
coordinates(krig_grid_small) <- ~ Var1 + Var2
#create kriging fit and apply interpolation to grid
krig_fit_small<-fields::Krig(ozone_df[1:2],ozone_df$Ozone)
ozone_krig_small<-raster::interpolate(grid_raster_small, krig_fit_small)
crs(ozone_krig_small) <-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
#plot output raster
plot(ozone_krig_small)

leaflet() %>% 
  addRasterImage(ozone_krig_small, project=T)%>% 
  addTiles() %>% 
  setView(lng = -111.941004,  40.610497, zoom = 10) %>% 
  addMiniMap()

Even though this plots with plot, when I try to add it to the leaflet I get Error in wkt(projfrom) : could not find function "wkt", which seems to be because I haven't set up the coordinates correctly for the raster image.


Solution

  • That error should go away if you update your packages (or at least raster, sp and rgdal) --- it is also good to use R >= 4.

    Here is a version with a few fixes that works for me

    library(raster)
    library(leaflet)
    library(fields)
    ozone_df<-data.frame("Longitude"=runif(200, -112.245075*10000,-111.455581*10000)/10000)
    ozone_df$Latitude<-runif(200, 40.063614*10000,40.827281*10000)/10000
    ozone_df$Ozone<-runif(200, 0,115)
    
    #create grid tick marks
    small_grid_x = seq(-111.455581,-112.245075,length.out=500)
    small_grid_y = seq(40.063614,40.827281,length.out=500)
    
    #create grid nodes
    krig_grid_small <- expand.grid(small_grid_x,small_grid_y)
    coordinates(krig_grid_small) <- ~ Var1 + Var2
    
    #create kriging fit and apply interpolation to grid
    krig_fit_small<-fields::Krig(ozone_df[1:2],ozone_df$Ozone)
    
    grid_raster_small = rasterFromXYZ(krig_grid_small)
    ozone_krig_small<-raster::interpolate(grid_raster_small, krig_fit_small)
    crs(ozone_krig_small) <-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
    #plot output raster
    plot(ozone_krig_small)
    
    leaflet() %>% 
      setView(lng = -111.941004,  40.610497, zoom = 10) %>% 
      addTiles() %>% 
      addRasterImage(ozone_krig_small, project=T)%>% 
      addMiniMap()