时间戳

首页 > 解决方案 > 如何设置插值栅格的投影以在 r 的传单图中绘制？

r - 如何设置插值栅格的投影以在 r 的传单图中绘制？

问题描述

我正在尝试将光栅图像添加到传单。虽然对我的插值并不完全满意，但我确实有一个通过正常plot函数绘制的图。我在下面显示的代码应该能够独立运行，因为我为输入数据提供了示例数据。

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()

即使这个情节带有情节，当我尝试将它添加到传单时，我得到了Error in wkt(projfrom) : could not find function "wkt"，这似乎是因为我没有为光栅图像正确设置坐标。

标签： rleafletraster

解决方案

如果你更新你的包（或至少 raster、sp 和 rgdal），这个错误应该会消失 --- 使用 R >= 4 也很好。

这是一个包含一些对我有用的修复的版本

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()

推荐阅读