微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

从脚本到 R

如何解决从脚本到 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 举报,一经查实,本站将立刻删除。