spatial - How to extract from a .nc file based on a shapefile in R? -


i have .nc file has global data , want extract data within boundary of .shp file. have tried several methods, still have problems.

the .nc file can downloaded @ https://www.dropbox.com/s/0ba6v4fnck8wjdm/spei03.nc?dl=0 , .shp file can downloaded @ https://www.dropbox.com/s/8wfgf8207dbh79r/gpr_000b11a_e.zip?dl=0

library(rgdal) library(ncdf4) library(raster)  shpfile<-readogr("gpr_000b11a_e.shp", layer="gpr_000b11a_e") g <- sptransform(shpfile, crs("+proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=nad83 +units=m +no_defs +ellps=grs80 +towgs84=0,0,0")) ncdata = raster(x="spei03.nc",varname="spei") proj4string(ncdata) = proj4string(g) mydata = rastertopoints(mask(x=ncdata,mask = g)) 

but cannot data. help.

to extract values, use extract() function. but, need select time step want use:

library(ncdf4) library(raster)  shpfile<-shapefile("gpr_000b11a_e.shp") # load shapefile  spei_nc=nc_open("spei03.nc") # open .nc data mtrx <- ncvar_get(spei_nc, varid="spei")[,,3] # extact value. in case, 3 time step number 3  spei_r <- raster(t(mtrx),xmn=-180, xmx=180, ymn=-90, ymx=90) # create raster data (it flipped) spei_r <- flip(spei_r,2) # correct raster projection(spei_r) <- '+proj=longlat +datum=wgs84 +ellps=wgs84 +towgs84=0,0,0 ' # add crs  plot(spei_r) # check result 

enter image description here

# reproject shapefile g <- sptransform(shpfile, spei_r@crs)  # extract values ecah feature spei_values <- extract(spei_r,g)  # extract mean each feature spei_mean <- extract(spei_r,g,fun=mean,na.rm=t)  spei_mean ##              [,1] ##  [1,] -0.37517873 ##  [2,]  0.24568238 ##  [3,]  0.12212451 ##  [4,] -1.44496705 ##  [5,]  0.57723304 ##  [6,]  0.08963697 ##  [7,]  0.07029936 ##  [8,]  0.01176584 ##  [9,]  0.35061048 ## [10,]  0.09197929 ## [11,]  0.41005611 ## [12,]  0.19130312 ## [13,] -0.69751740 

Comments

Popular posts from this blog

ZeroMQ on Windows, with Qt Creator -

unity3d - Unity SceneManager.LoadScene quits application -

python - Error while using APScheduler: 'NoneType' object has no attribute 'now' -