usmap.Rmd
require(raster)
require(sp)
require(ggplot2)
require(PEcAn.visualization)
data(yielddf, package = "PEcAn.visualization")
#pecan.worldmap(yielddf, outfile=file.path(tempdir(), 'foo.png'))
spdf <- SpatialPointsDataFrame( data.frame( x = testdf$y , y = testdf$x ) , data = data.frame( z = testdf$z ) )
# Plotting the points reveals the unevenly spaced nature of the points
spplot(spdf)
files <- dir("~/dev/bety/local/modelout", pattern="grid.csv", full.names=TRUE)
yieldfiles <- files[!grepl("evapotranspiration", files)]
etfiles <- files[grepl("evapotranspiration", files)]
for(file in yieldfiles){
df.in <- read.csv(file)
outfile <- gsub("csv", "png", file)
#pecan.worldmap(df.in, outfile=outfile)
}
for(file in etfiles){
df.in <- read.csv(file, skip = 1)
outfile <- gsub("csv", "png", file)
#pecan.worldmap(df.in, outfile=outfile)
}
# Make an evenly spaced raster, the same extent as original data
e <- extent( spdf )
# Determine ratio between x and y dimensions
ratio <- ( e@xmax - e@xmin ) / ( e@ymax - e@ymin )
# Create template raster to sample to
r <- raster( nrows = 56 , ncols = floor( 56 * ratio ) , ext = extent(spdf) )
rf <- rasterize( spdf , r , field = "" z, fun = mean )
# Attributes of our new raster (# cells quite close to original data)
rf
# We can then plot this using `geom_tile()` or `geom_raster()`
rdf <- data.frame( rasterToPoints( rf ) )
ggplot( NULL ) + geom_raster( data = rdf , aes( x , y , fill = layer ) )
# from http://gis.stackexchange.com/a/20052/3218
require(rgdal)
proj4string(spdf) <- CRS("+init=epsg:4326")
pts <- spTransform(spdf ,CRS("+proj=longlat +datum=WGS84"))
gridded(spdf) <- TRUE
p2g <- points2grid(spdf, tolerance = 0.95)
r <- raster(p2g)
plot(p2g)
this was run in data from Willow manuscript. See redmine issue 2387
files <- dir(".", pattern = ".txt")
library(data.table)
lat <- fread(files[1])[,list(lat = V2)]$lat
lon <- fread(files[1])[,list(lon = V3)]$lon
system.time(a <- fread(files[1]))
#yield <- data.table(date = vector(), lat = vector(), lon = vector(), yield = vector())
yieldarray <- array(NA, dim = c(56, 135, 0), dimnames = c('lat', 'lon', 'year'))
for(file in files){
newyield <- fread(file)[,list(lat = V2, lon = V3, yield = V24)][!is.na(yield)]
newyield.regrid <- regrid(newyield)
yieldarray <- abind(yieldarray, newyield.regrid)
}
y <- list()
for(year in 1979:2010){
file <- files[grepl(year, files)]
y[[as.character(year)]] <- fread(file)[,list(yield = V24)]
}
a <- do.call("cbind", y)
setnames(a, as.character(1979:2010))
b <- a[,list(yield = rowMeans(.SD))]$yield
willow <- data.table(lat = lat, lon = lon, yield = b)[!is.na(yield)]