0

I have a script which reads a .dat file containing co-ordinates, rainfall and crop area data as below.

  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 

The code below calculates the emissions for Wheat through a series of models before extracting the data and producing a raster file and plotting it in ggplot. Having read these related queries I'm getting rather confused as to whether I need to make a package or am missing something very basic as to how to repeat the code through each crop type.

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

I'm trying to turn as much of the above into functions as I can but they're all rather more complicated than the examples I find on here and in the ?help files. Many thanks in advance for any help/suggestions.

B.Wel
  • 86
  • 10

1 Answers1

1

You need to put your code in a function that takes the data file name as argument. Then you can call your function with the relevant data file name. Something like:

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

If the model parameters for different data, you need to add the vector of parameter values as argument to the function.

citronrose
  • 409
  • 3
  • 9
  • Cool, thanks! If I subset the data into one for each species (Wheat, SBarley, etc.) if I use `crm[[3]]` in place of `crm$Wheat` will that read the third column in each subset dataframe? – B.Wel Aug 19 '16 at 07:45
  • Realized the easiest fix to get the function to work was to simplify the original data frame to just five columns using the `melt` function from `reshape2`. – B.Wel Aug 22 '16 at 08:11