如何解决从脚本到 R
我正在处理大学关于 Rstudio 的练习
我在 R 中编写了一个脚本,该脚本基于 csv 文件执行克里金插值,该脚本中不包含函数,我使用了几个库中的一些函数。
脚本运行良好,最后我得到了一个带有克里金表面的光栅文件和来自该光栅的图形。
但是现在我必须从该脚本创建一个包,我按照互联网上的一些教程进行了尝试,但它不起作用我认为是因为我的脚本不是函数
我应该怎么做才能从我的脚本创建一个包?
library(sf)
library(sp)
library(raster)
library(rgdal)
library(gstat)
library(tmap)
#configuring folder for the result
setwd("C:")
dir.create("R_Exercise")
setwd("R_Exercise")
getwd()
#reading csv from google drive
id <- "1YOCFoHkKqkc69YfvdV-lrUThTP_32LHn" # google file ID
read.csv(sprintf("https://docs.google.com/uc?id=%s&export=download",id))
silver_data = read.csv(sprintf("https://docs.google.com/uc?id=%s&export=download",id))
#checking imported data
str(silver_data)
#data exploration
silver_data$silver
hist(silver_data$silver,breaks = 10)
summary(silver_data$silver)
#Creating Spatial sf object
silver_2 <- st_as_sf(silver_data,coords = c("Este","norte"))
#checking the spatial object
str(silver_2)
#asigning crs from epsg codes: 9155 SIRGAS-Chile 2016 / UTM zone 19S
silver_2 <- st_set_crs(silver_2,9155)
st_crs(silver_2)
#ploting silver data
plot(silver_2["silver"],asp = 1,cex = 4*silver_2$silver/max(silver_2$silver),main= "Silver Data")
#Starting to solve the problem find the kriging surface for silver grade
#empirical variogram
varem <- variogram(silver ~ 1,silver_2,cutoff = 5000,width = 200)
plot(varem,plot.numbers = T,asp=1)
#experimental variogram
varex <- vgm(psill = 350,model = "Sph",range = 2000,nugget = 0)
plot(varem,pl = T,model = varex)
#automatic adjustment
varaut <- fit.variogram(varem,varex)
plot(varem,model = varaut)
#interpolation
#creating empty grid
emptyr <- raster(as(silver_2,"Spatial"),ncols = 200,nrows = 80)
#cheking empty raster
extent(emptyr)
#resize spatial resolution
res(emptyr)<-c(100,100)
Emptygrid <- rasterToPoints(emptyr,spatial = TRUE)
gridded(Emptygrid) <- TRUE
Emptygrid <- as(Emptygrid,"SpatialPixels")
Emptygrid <- sf::st_as_sf(Emptygrid)
#obtaining kriging
Inter <- krige(silver~1,locations=silver_2,newdata=Emptygrid,model=varaut)
Inter
#ploting some verifications
Inter_sp <- as(Inter,Class = "Spatial")
print(spplot(Inter_sp,"var1.pred",asp=1,main="Prediction for Silver Grade,gr/Ton"))
print(spplot(Inter_sp,"var1.var",col.regions=cm.colors(64),main="Prediction variance"))
Inter_sp$coeffvar = sqrt(Inter_sp$var1.var)/Inter_sp$var1.pred
print(spplot(Inter_sp,"coeffvar",main="Variation Coefficient"))
#creating final file
r <-raster(Inter_sp)
res(r)<-c(100,100)
img <-rasterize (x = Inter_sp,y = r,field = Inter_sp$var1.pred )
writeraster(img,'output.tif',overwrite = TRUE)
result <-raster("output.tif")
plot(result)
#Final Graphic
tm_shape(result) + tm_raster(n=10,palette="RdBu",auto.palette.mapping=FALSE,midpoint = NA,title="Kriging Prediction \nSilver grade \n (gr/Ton)") + tm_shape(silver_2) + tm_dots(size=0.2) + tm_legend(legend.outside=TRUE)
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。