[R-sig-Geo] Translate a net-cdf file with different years to polygon regions.

Stachelek, Joseph jstachel at sfwmd.gov
Tue Apr 19 14:42:51 CEST 2016

Since your nc file contains multiple layers, you will want to use `raster::stack()` rather than `raster::raster()`.

-----Original Message-----
From: Stachelek, Joseph
Sent: Tuesday, April 19, 2016 8:23 AM
To: 'Janka VANSCHOENWINKEL' <janka.vanschoenwinkel at uhasselt.be>; r-sig-geo at r-project.org
Subject: RE: [R-sig-Geo] Translate a net-cdf file with different years to polygon regions.

Hi Janka,

I think you can simplify your code a lot by opening your nc file directly using the `raster` package rather than messing with `nc_open` calls.

ncin <- raster::raster(paste(ncname, ".nc", sep = ""))

Then you might use `raster::extract()` to pull the values associated with your polygons. Also, I would recommend posting a link to a gist (https://gist.github.com/) rather than pasting such a long script into your email.


-----Original Message-----
From: R-sig-Geo [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of Janka VANSCHOENWINKEL
Sent: Tuesday, April 19, 2016 4:45 AM
To: r-sig-geo at r-project.org
Subject: [R-sig-Geo] Translate a net-cdf file with different years to polygon regions.


I have a net-cdf file which has raster of 0.5 on 0.5 degrees. You can
easily download it here
search for *cru_ts3.23.2001.2010.pet.dat.nc.gz *(it is also downloadable as
a dat.file if this is more handy to work with)
 (Or simply download the net-cdf file directly through:

I opened the file in ArcMap as well and found that the coordinate system
used is: GCS_WGS_1984. The net-cdf file contains monthly data from 2001-2010

Download from this website
following ZIP file: *NUTS_2006_60M_SH.zip* and save all the *NUTS_RG_60M_2006
files *in a folder where you can easily find the back. In what follows I
refer to this as "NUTS3".


- I want to add the information from the raster files to the NUTS3
shapefile (which are polygons and no rasters) in order to obtain a table
per nuts3 region for each monthtly variable.


The file appears to be very difficult to work with in ArcGis. Also, I will
have to repeat this process a number of times for different variables. So I
would like to have one code in R that I can run.

I have a number of separate codes to do the following:
- translate the net-cdf file to a cvs file with longitude and latitude as
identifiers (see below under 1)
- translate a cvs file with accompanying empty raster file to nuts3
regions. (this is a code that I have used before when I had a cvs file and
a raster). (see below under 2).

However, I don't have a raster file now. Well, technically I do (since the
net-ncf file is a raster) but I don't know how to use it in this format.

Can somebody help me to link the codes below or suggest a different code to
obtain what I want?

Thanks a lot!


*1) With the following code, I can make a cvs file and extract all the data
in table format.*


ncfname <- paste(ncname, ".nc", sep = "")
dname <- "pet"
ncin <- nc_open(ncfname)

lon <- ncvar_get(ncin, "lon")
nlon <- dim(lon)

lat <- ncvar_get(ncin, "lat", verbose = F)
nlat <- dim(lat)

print(c(nlon, nlat))

t <- ncvar_get(ncin, "time")
tunits <- ncatt_get(ncin, "time", "units")
nt <- dim(t)

tmp.array <- ncvar_get(ncin, dname)
dlname <- ncatt_get(ncin, dname, "long_name")
dunits <- ncatt_get(ncin, dname, "units")
fillvalue <- ncatt_get(ncin, dname, "_FillValue")

title <- ncatt_get(ncin, 0, "title")
institution <- ncatt_get(ncin, 0, "institution")
datasource <- ncatt_get(ncin, 0, "source")
references <- ncatt_get(ncin, 0, "references")
history <- ncatt_get(ncin, 0, "history")
Conventions <- ncatt_get(ncin, 0, "Conventions")


# split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth = as.integer(unlist(tdstr)[2])
tday = as.integer(unlist(tdstr)[3])
tyear = as.integer(unlist(tdstr)[1])
chron(t, origin = c(tmonth, tday, tyear))

tmp.array[tmp.array == fillvalue$value] <- NA

length(na.omit(as.vector(tmp.array[, , 1])))

m <- 1
tmp.slice <- tmp.array[, , m]
image(lon, lat, tmp.slice, col = rev(brewer.pal(10, "RdBu")))

grid <- expand.grid(lon = lon, lat = lat)
cutpts <- c(-50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50)
levelplot(tmp.slice ~ lon * lat, data = grid, at = cutpts, cuts = 11,
pretty = T,
          col.regions = (rev(brewer.pal(10, "RdBu"))))

lonlat <- expand.grid(lon, lat)
tmp.vec <- as.vector(tmp.slice)

tmp.df01 <- data.frame(cbind(lonlat, tmp.vec))
names(tmp.df01) <- c("lon", "lat", paste(dname, as.character(m), sep = "_"))
head(na.omit(tmp.df01), 20)

csvfile <- "cru_tmp_1.csv"
write.table(na.omit(tmp.df01), csvfile, row.names = FALSE, sep = ",")

tmp.vec.long <- as.vector(tmp.array)

tmp.mat <- matrix(tmp.vec.long, nrow = nlon * nlat, ncol = nt)


lonlat <- expand.grid(lon, lat)
tmp.df02 <- data.frame(cbind(lonlat, tmp.mat))

names(tmp.df02) <- c("lon","lat","pet_jan_2001",

options(width = 110)
head(na.omit(tmp.df02, 20))


csvfile <- "cru_tmp_2.csv"
write.table(na.omit(tmp.df02), csvfile, row.names = FALSE, sep = ",")

*2) translate a cvs-file with accompanying raster file to polygon regions.*

The  "filename.txt" file should contain the variables: lon, latitude, and
all the monthly_yearly variables extracted from point 1 above.

The grid shapefile (*grid_025dd.shp*) can be found through the following
link but it is only an example and not the correct grid for the problem
above :

# upload data
mydata<-read.table("filename.txt", header=TRUE,sep=",",dec=".")

# upload empty raster
# 40 seconds
grid <- readOGR(".", layer = "grid_025dd")

# concatenate data in R
# 2 seconds
mydata$lonlat<-do.call(paste, c(mydata[c("lon", "lat")], sep=""))
grid at data$lonlat<-do.call(paste, c(grid at data[c("LONGITUDE", "LATITUDE")],

# use common variable lonlat to merge data in raster

###### prepare shapefile #####
library(rgdal)        ## Load geographic info
library(maps)         ## Projections
library(maptools)     ## Data management
#library(sm)           ## Data management
library(spdep)        ## Spatial autocorrelation
library(gstat)        ## Geostatistics
library(splancs)      ## Kernel Density
library(spatstat)     ## Geostatistics
library(pgirmess)     ## Spatial autocorrelation
library(RColorBrewer) ## Visualization
library(classInt)     ## Class intervals
library(spgwr)        ## GWR

# Match polygons with data
idx <- match(grid$lonlat, mydata$lonlat)
# Places without information
idxNA <- which(is.na(idx))
# Information to be added to the SpatialPolygons object
dat2add <- mydata[idx, ]
# spCbind uses row names to match polygons with data
# First, extract polygon IDs
IDs <- sapply(grid at polygons, function(x)x at ID)
# and join with the SpatialPolygons
row.names(dat2add) <- IDs
datPols <- spCbind(grid, dat2add)
# Drop those places without information
datPols <- datPols[-idxNA, ]
# write new shapefile
# 7 seconds
writeOGR(datPols, dsn = ".", layer ='sm2000eu28', driver = 'ESRI Shapefile')
# read new shapefile
# 51 seconds
data <- readOGR(".", layer="sm2000eu28")

# intersect nuts with grid #

nuts <- readOGR(".", layer = "NUTS_RG_60M_2006")

proj4string(data) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")

grid <- data
grid at data$lonlat <- NULL
grid at data$lonlat_1 <- NULL
grid at data$ID <- NULL
grid at data$lat <- NULL
grid at data$lon <- NULL
grid at data$ELEVATION <- NULL
grid at data$DAYS_RAIN_ <- NULL

# First find out which grid cells intersect your NUTS polygons
grid_nuts <- gIntersects(grid,nuts,byid = TRUE)

# use the apply() function to calculate the mean, min, and max of your
# The loop makes

for(i in names(grid at data)){
  nuts at data[[paste(i, 'average_value', sep="_")]] <-
apply(grid_nuts,1,function(x) mean(grid at data[[i]][x]))
  nuts at data[[paste(i, 'min_value', sep="_")]] <-
apply(grid_nuts,1,function(x) min(grid at data[[i]][x]))
  nuts at data[[paste(i, 'max_value', sep="_")]] <-
apply(grid_nuts,1,function(x) max(grid at data[[i]][x]))

write.table(nuts at data, "nuts_sm2000eu28_unweighted.txt", sep="\t")

        [[alternative HTML version deleted]]

R-sig-Geo mailing list
R-sig-Geo at r-project.org

We value your opinion. Please take a few minutes to shar...{{dropped:5}}

More information about the R-sig-Geo mailing list