Extract values at XY position from NetCDF file

ncdf_extract_values.R

library(raster)
library(ncdf4)
library(maptools)

# netcdf file from ftp://aspen.atmos.albany.edu/pdsi/cmip5/scPDSIpm/pdsisc.monthly.maps.1900-2099.r2.5x2.5.EnsAvg14Models.TP2.ipe=2.nc
# shp file from http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/ne_50m_airports.zip

# load netcdf file as a raster brick
f <- "D:/bacasable/pdsisc.monthly.maps.1900-2099.r2.5x2.5.EnsAvg14Models.TP2.ipe=2.nc"
rbr <- brick(f, varname="pdsisc")
# display title + nb layers
rbr@title
rbr@data@nlayers
# get layers name
lyr_names <- rbr@data@names

# plot layers 1, 6, 1200, 1201
plot(raster(rbr,layer=1), main=lyr_names[1])
plot(raster(rbr,layer=6), main=lyr_names[6])
plot(raster(rbr,layer=1200), main=lyr_names[1200])
plot(raster(rbr,layer=1201), main=lyr_names[1201])

# read point shapefile : airport
pts <- readShapePoints("D:/GIS_DATA/NaturalEarth/50m_cultural/ne_50m_airports.shp")
# get xy coords as a matrix + get iata_code field from attribute table + display first values
mat_coords <- pts@coords
iata_code <- pts@data$iata_code
head(mat_coords)

# extract pdsisc values (2400) at points positions (281)
# we get a 281 rows x 2400 cols matrix
pdsisc_extract <- extract(rbr, mat_coords)
# transpose matrix => 2400 rows x 281 cols + set iata_code as cols name
pdsisc_extract <- t(pdsisc_extract)
dimnames(pdsisc_extract)[[2]] <- iata_code

# bind layers name and values + save data in a tab separated file
df_out <- cbind(lyr_names, pdsisc_extract)
write.table(df_out, file="D:/bacasable/pdsisc_extract.txt", sep="\t", dec=",", row.names=F)