require(raster)
require(ncdf4)
require(sp)
require(rgdal)
setwd("C:/Users/jvithana/Documents/MobaXterm/home")
stdy<- readOGR(dsn = "C:/Users/jvithana/Documents/MobaXterm/home", layer = "MIAMI_dade-final")
stdy<- spTransform(stdy, CRS("+init=epsg:4326"))
rs<- stack("C:/Users/jvithana/Documents/MobaXterm/home/wrfout_d03_2005-08-24_00100100.nc", varname = "RAINNC")
strsplit(rs@layers[[3]]@file@name, "_")[[1]][3:5]
Extract the date from the file name
dt<- paste0(strsplit(rs@layers[[3]]@file@name, "_")[[1]][3:5], collapse = "1")
vars<- c("RAINNC")
outDF<- data.frame()
for(i in 1:length(vars)){
r<- stack("C:/Users/jvithana/Documents/MobaXterm/home/wrfout_d03_2005-08-24_00100100.nc", varname = vars[i])
Set the extent and crs
extent(r)<- c(-82.83185, -77.34158, 23.34295, 28.16767)
proj4string(r)<- "+init=epsg:4326"
Reduce raster to the study area
r<- crop(r, stdy)
outDF<- rbind(outDF, data.frame(Variable = vars[i], data.frame(r[c(39:30)])))
}
Coerce to a POSIXct
dt<- as.POSIXct(dt, format = "%Y-%m-%d", tz = "GMT")
Number of 6 hour increments
hr<- (seq(from = 0, to = (nlayers(r) - 1) * 1, by = 1)) * 60 * 60
Calculate the actual dates
dt<- dt + hr
dt<- as.character(dt)
Rename the columns
names(outDF)[2:ncol(outDF)]<- dt
#set the working directory
#setwd("/glade/work/jvithana/WRF60/out")
write.table(outDF, "wrfout_d03_2005-08-24_00100100.txt", sep = "\t", row.names=FALSE )
Crop and mask
#plot(rs)
plot(r[[90]])
#plot(stdy, add = TRUE)
Question,
my cropped grid consist with a 3D array 3039111 and it reads only 10 and I cannot figure out why it is reading only 10 pixels as I did not mention anything about 10 pixels.