以下のように、座標、降雨量、作付面積のデータを含む.datファイルを読み取るスクリプトがあります。
East North rain Wheat Sbarley Potato OSR fMaize Total LCA
10000 510000 1498.73 0.0021 0.5 0.0022 0 0.0056 0.01 0.01
10000 510000 1498.73 0.0021 0.034 0.0022 0 0.0056 0.01 0.01
10000 510000 1498.73 0.0021 0.001 0.0022 0 0.0056 0.01 0.01
10000 515000 1518.51 0.0000 0.12 0.0000 0 0.0000 0.00 0.00
10000 515000 1518.51 0.0000 0.0078 0.0125 0 0.0000 0.00 0.00
10000 515000 1518.51 0.0000 0 0.0000 0 0.03 0.00 0.00
以下のコードは、データを抽出してラスターファイルを作成し、それをggplotにプロットする前に、一連のモデルを通じて小麦の排出量を計算します。これらの関連するクエリを読んだ 後、パッケージを作成する必要があるのか、各クロップタイプでコードを繰り返す方法について非常に基本的なものが不足しているのかについて、かなり混乱しています。
library(doBy)
library(ggplot2)
library(plyr)
library(raster)
library(RColorBrewer)
library(rgdal)
library(scales)
library(sp)
### Read in the data
crmet <- read.csv("data.dat")
# Remove NA values
crm <- crmet[ ! crmet$rain %in% -119988, ]
crm <- crm[ ! crm$Wheat %in% -9999, ]
### Set model parameters
a <- 0.1474
b <- 0.0005232
g <- -0.00001518
d <- 0.000003662
N <- 182
### Models
crm$logN2O <- a+(b*crm$rain)+(g*N)+(d*crm$rain*N)
crm$eN2O <- exp(crm$logN2O)
crm$whN2O <- crm$eN2O*crm$Wheat
### Prepare data for conversion to raster
crmet.ras <- crm
crmet.ras <- rename(crmet.ras, c("East"="x","North"="y"))
#### Make wheat emissions raster
wn <- crmet.ras[,c(1,2,13)]
spg <- wn
# Set the Eastings and Northings to coordinates
coordinates(spg) <- ~ x + y
# coerce to SpatialPixelsDataFrame
gridded(spg) <- TRUE
# coerce to raster
rasterDF <- raster(spg)
# Add projection to it - in this case OSBG36
proj4string(rasterDF) <- CRS("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs")
rasterDF
writeRaster(rasterDF, 'wn.tif', overwrite=T)
### Plotting the raster:
whplot <-ggplot(wn, aes(x = x, y = y))+
geom_tile(aes(fill = whN2O))+
theme_minimal()+
theme(plot.title = element_text(size=20, face="bold"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.title = element_text(size=16, face="bold"))+
scale_fill_gradient(name=expression(paste("", N[2], "O ", Ha^-1, sep="")))+
xlab ("")+
ylab ("")+
labs (title="Nitrous oxide emissions\nfrom Wheat")
whplot
私は上記の多くをできる限り関数に変えようとしていますが、それらはすべて、ここや?help
ファイルにある例よりもかなり複雑です。助け/提案を事前に感謝します。
データファイル名を引数として取る関数にコードを配置する必要があります。次に、関連するデータファイル名を使用して関数を呼び出すことができます。何かのようなもの:
library(doBy)
library(ggplot2)
library(plyr)
library(raster)
library(RColorBrewer)
library(rgdal)
library(scales)
library(sp)
#defining the function
my.neat.function <- function(datafname){
### Read in the data
crmet <- read.csv(datafname)
# Remove NA values
crm <- crmet[ ! crmet$rain %in% -119988, ]
crm <- crm[ ! crm$Wheat %in% -9999, ]
### Set model parameters
a <- 0.1474
b <- 0.0005232
g <- -0.00001518
d <- 0.000003662
N <- 182
### Models
crm$logN2O <- a+(b*crm$rain)+(g*N)+(d*crm$rain*N)
crm$eN2O <- exp(crm$logN2O)
crm$whN2O <- crm$eN2O*crm$Wheat
### Prepare data for conversion to raster
crmet.ras <- crm
crmet.ras <- rename(crmet.ras, c("East"="x","North"="y"))
#### Make wheat emissions raster
wn <- crmet.ras[,c(1,2,13)]
spg <- wn
# Set the Eastings and Northings to coordinates
coordinates(spg) <- ~ x + y
# coerce to SpatialPixelsDataFrame
gridded(spg) <- TRUE
# coerce to raster
rasterDF <- raster(spg)
# Add projection to it - in this case OSBG36
proj4string(rasterDF) <- CRS("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs")
rasterDF
writeRaster(rasterDF, 'wn.tif', overwrite=T)
### Plotting the raster:
whplot <-ggplot(wn, aes(x = x, y = y))+
geom_tile(aes(fill = whN2O))+
theme_minimal()+
theme(plot.title = element_text(size=20, face="bold"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.title = element_text(size=16, face="bold"))+
scale_fill_gradient(name=expression(paste("", N[2], "O ", Ha^-1, sep="")))+
xlab ("")+
ylab ("")+
labs (title="Nitrous oxide emissions\nfrom Wheat")
whplot
} #end of function definition
my.neat.function("data.dat") #first call to function
my.neat.function("otherdata.dat")#same thing with another dataset
異なるデータのモデルパラメーターの場合、パラメーター値のベクトルを引数として関数に追加する必要があります。
この記事はインターネットから収集されたものであり、転載の際にはソースを示してください。
侵害の場合は、連絡してください[email protected]
コメントを追加