1

I have managed to open and intersect a gridded .nc file with spatial polygons for river basins. Unfortunately, I cannot find which variable indicates the month at when the values were measured. The file should contain monthly irrigated area data.

Here is my code to open the file

irrig_path <- "/Volumes/Transcend/Uni/ETH/International Political Economy/data/Huang Water Use/irrigation water use/"
irrig_name <- "cons_irr_h08.nc"
irrig_file <- paste(irrig_path, irrig_name, sep="")
irrig_r <- raster(irrig_file, varname="cons_irr")
irrig_r

Here is a link to the file

https://drive.google.com/file/d/1Ry68E_M3e0twpBn8V8NdwSsHKkERE3AR/view?usp=sharing

and a link to the article that published the data. I could not find any indication where the date is stored.

https://www.hydrol-earth-syst-sci.net/22/2117/2018/hess-22-2117-2018.pdf 

This is an image of the content of the raster data once loaded into R

enter image description here

0

Looking at the netcdf file with ncdump -h cons_irr_h08.nc indicates the dimensions of the file are 480 months by 64028 grid_num, with a level of abstraction from grid_num into lat, lon. The data is structured as a set of 64028 individual time series. I don't think this will easily parse into the raster package.

$ ncdump -h cons_irr_h08.nc 
netcdf cons_irr_h08 {
dimensions:
    grid_num = 64028 ;
    month = 480 ;
variables:
    double cons_irr(month, grid_num) ;
        cons_irr:units = "mm/month" ;
        cons_irr:description = "1. gridded irrigation consumption results: 64028 rows, 480 Months;2.the results is based on the outputs of H08 forcing by WFDEI" ;
    double lon(grid_num) ;
        lon:units = "degree" ;
    double lat(grid_num) ;
        lat:units = "degree" ;
    double month(month) ;
        month:units = "months since 1971-1" ;
}

Maybe use library(ncdf4) or library(RNetCDF) to deal with file more directly and consider it as spatial point data.

ETA:

Looking at the lat and lon variables with library(ncdf4) shows it is pretty sparse and you need some sort of interpolation tool to re-grid the data to a raster.

library(ncdf4)
ncin <- nc_open("cons_irr_h08.nc")

lon= ncvar_get(ncin,'lon')
lat= ncvar_get(ncin,'lat')

head(cbind(lat,lon))
       lat   lon
[1,] 32.75 60.75
[2,] 33.25 60.75
[3,] 33.75 60.75
[4,] 34.25 60.75
[5,] 29.75 61.25
[6,] 30.25 61.25

plot(lon,lat,type='p',pch='.')

plot of datapoint locations

The "remap the entire field" section from https://rpubs.com/markpayne/132500 or https://gis.stackexchange.com/a/35322/10229 or How to make RASTER from irregular point data without interpolation might be of some help with this data.

Back to the time variable question, you can extract it directly with the ncdf4 package:

library(ncdf4)
ncin <- nc_open("cons_irr_h08.nc")    
tvar= ncvar_get(ncin,'month')

... and then you could turn the numerics into dates with this trick adapted from https://gis.stackexchange.com/a/35322/10229 :

add.months= function(date,n) seq(date, by = paste (n, "months"), length = 2)[2]
times = as.Date(sapply(tvar-1,add.months,date=as.Date('1971-01-15')),origin='1970-01-01')

head(times)

[1] "1971-01-15" "1971-02-15" "1971-03-15" "1971-04-15" "1971-05-15" "1971-06-15"
  • Looking a little bit at the lat,lon variables, they do not seem to be in any sort of raster order. I don't know the dataset, but to make a raster from random points, you have to do some sort of interpolation. – Dave X Nov 26 '18 at 16:13
  • Your help is highly appreciated Dave X! If I now retrieve the variable cons_irr cons_irr <- ncvar_get(ncin, 'cons_irr') Couldn't I now just simply make a spatialPointDataframe form the extracted lon, lat, times and cons_irr? If so, how could I do that? – S Front Nov 26 '18 at 21:04
  • I found this wesite here that shows how to do it but it doesn't show how to do it for several observations over time rspatial.org/spatial/rst/3-vectordata.html Can I assume that the first entry in cons_irr can be associated with the first two entries of cbind(lon, lat)? – S Front Nov 26 '18 at 21:08
  • I don't have the dataset handy anymore to test, but I think you can try n=1;xx=nc_varget(ncin,'cons_irr',start=c(n,1),count=c(1,-1)) to get the full time series associated with the first (n=1) point. Or n=1;xx=nc_varget(ncin,'cons_irr',start=c(1,n),count=c(-1,1)) to get all the points at times step n=1. – Dave X Nov 27 '18 at 4:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.