I'm fairly new to using R for GIS purposes. I have a netcdf file containing several variables with multiple dimensions (x,y,z,value and time). I am trying to turn this into a raster brick. The data is quite large so I need to pull data from a specified time window and z(depth). This has not been a problem and extract an array with the appropriate dimensions using the below code .
library(ncdf4)
library(raster)
t <- ncvar_get(nc, "model_time")
t1<-ncvar_get(nc,"model_time_step")
tIdx<-t[t> 20120512 & t < 20120728]
tIdx2<-which(t> 20120512 & t < 20120728)
# Depths profiles < 6 meters
dIdx<-which(nc$dim$depthu$vals <6)
# ncdf dimension lengths
T3 <- nc$var[[7]]
varsize <- T3$varsize
# Define the data (depths,time,etc.) you wish to extract from the ncdf
start <- c(x = 1, y= 1,depthu=1, time_counter = min(tIdx2))
count <- c(x = max(varsize[1]), y = max(varsize[2]),depthu=1, time_counter =
max(tIdx2)-min(tIdx2)+1)
# order of the dimensions
dim.order <- sapply(nc$var$votemper$dim, function(x) x$name)
temp<-ncvar_get(nc,"votemper",start=start[dim.order],count=count[dim.order])
nc$var$votemper
An example of my data (dropping the depth/z and the time dimensions)
temp<-structure(c(0,0,0,0,0,0,0,15.7088003158569,15.3642873764038,14.9720048904419,,15.9209365844727,14.9940872192383,15.0184164047241,15.0260219573975, 0,15.7754755020142, 15.424690246582, 15.6697931289673,15.6437339782715, 0,15.6151847839355, 15.5979156494141, 15.6487197875977,15.432520866394), .Dim = c(x = 5L, y = 5L))
The latitudes and longitudes extracted from the ncdf are irregularly spaced and two dimensions each (i.e. An irregular spaced lat and lon for each cell)
lon<-structure(c(-71.2870483398438,-71.2038040161133,-71.1205596923828,-71.0373153686523, -70.9540710449219, -71.2887954711914, -71.2055587768555,-71.122314453125, -71.0390701293945,-70.9558258056641,-71.2905654907227,-71.2073211669922,-71.1240844726562,-71.0408401489258,-70.9576034545898,-71.292350769043,-71.209114074707, -71.1258773803711, -71.0426330566406,-70.9593963623047, -71.2941513061523, -71.2109222412109, -71.127685546875,-71.0444488525391, -70.9612045288086), .Dim = c(5L, 5L))
lat<-structure(c(38.5276718139648, 38.529125213623, 38.5305824279785,38.532039642334, 38.5334968566895, 38.5886116027832, 38.5900802612305,38.591552734375, 38.5930252075195, 38.5944976806641, 38.6494789123535,38.6509628295898, 38.6524467468262, 38.6539344787598, 38.6554222106934,38.7102699279785, 38.7117652893066, 38.713264465332, 38.7147674560547,38.7162704467773, 38.7709808349609, 38.7724952697754, 38.7740097045898,38.7755241394043, 38.777042388916), .Dim = c(5L, 5L))
Typically I would generate a raster brick from this data using
Temp_brick <- brick(temp, xmn=min(lat), xmx=max(lat), ymn=min(lon), ymx=max(lon),transpose=T)
Temp_brick<-t(flip(Temp_brick,1))
This, however does not account for the irregular spacing and raster cell values are located in the wrong position (lon,lat). I have searched across stack overflow and other gis help sources and I can't find a similar problem with a solution or I'm not asking the right question. I'm not particularly sure how to go about this. Not sure whether this should be dealt with when extracting the data from the netcdf or if it should be dealt with after the raster brick has been created without defined extent. I have tried to find a way to define the lon lats for the raster without any luck. Tried converting lon,lat and value to 3 column dataframe and then use the raster::rasterFromXYZ function. This won't work quick enough for the size of the data I'm dealing with, which in reality is 197(x)*234(y)*2(z)*900(time)*5(variables)*12(years(separate netcdf files).
Any help is greatly appreciated
an option with akima to first interp the data to a regular grid and then turn it into a raster:
# define the regular lon lat or just pass the nx, ny param to interp functions
lonlat_reg <- expand.grid(lon = seq(min(lon), max(lon), length.out = 5),
lat = seq(min(lat), max(lat), length.out = 5))
# interp irregular data to a regular grid
# both solution return the same results because
# i've define the regular grid as akima default
test <- interp(x = as.vector(lon), y = as.vector(lat), z = as.vector(temp),
xo = unique(lonlat_reg[,"lon"]), yo = unique(lonlat_reg[,"lat"]),
duplicate = "error", linear = FALSE, extrap = FALSE)
test <- interp(x = as.vector(lon), y = as.vector(lat), z = as.vector(temp),
nx = 5, ny = 5, linear = FALSE, extrap = FALSE)
# turn into a raster
test_ras <- raster(test)
Check the arguments of the function to choose the interpolation performed etc and be careful if you use extrapolation!
I've seen also that method
Cheers