WRF Outputs in Netcdf format,

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.

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.