Translate a net-cdf file with different years to polygon regions.
On Fri, 29 Apr 2016 at 22:23 Janka VANSCHOENWINKEL <
janka.vanschoenwinkel at uhasselt.be> wrote:
Hi Joseph, Thank you very much for the hint. If I may ask you, could you please help me a little bit further. I am not familiar with these packages and I do not see how to overlap the point data or the netcdf data with the nuts3 polygons. Basically, I have too options: 1) use the cvs-points (with lon, lat, and the data variable) to see which points are located in which nuts3 regions and then save these results. (but then my previous codes still applies) 2) use the original netcdf file and overlay it with the nuts3 polygon shapefile (this was your hint, but I don't know how to start with this).
Something like this will do it:
library(raster)
require(rgdal)
require(ncdf4)
st <- stack("thencfile.nc") ## might need varname = "sst" or similar
poly <- shapefile("theshpfile.shp")
extract(st, poly, fun = "mean") ## and maybe na.rm = TRUE
If you don't want the mean, you can just leave out the fun argument and
you'll get all the pixel values for every time step in a big list, which
may not be obviously helpful but it's all useable with generic R.
I can't quite bring myself to get your files and try it, so please try and
report details.
Cheers, Mike.
Thank you very much for your answer! 2016-04-19 14:42 GMT+02:00 Stachelek, Joseph <jstachel at sfwmd.gov>:
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. Joe -----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. *DATA* I have a net-cdf file which has raster of 0.5 on 0.5 degrees. You can easily download it here <
and 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: cru_ts3.23.2001.2010.pet.dat.nc.gz <
). 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 <
the 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". *WHAT I WANT* - 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. *WHERE I AM STUCK* 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! Janka *1) With the following code, I can make a cvs file and extract all the
data
in table format.*
library(fields)
library(chron)
library(ncdf4)
ncname<-"cru_ts3.22.2001.2010.pet.dat"
ncname<-"cru_ts3.23.1991.2000.pet.dat"
ncfname <- paste(ncname, ".nc", sep = "")
dname <- "pet"
ncin <- nc_open(ncfname)
print(ncin)
lon <- ncvar_get(ncin, "lon")
nlon <- dim(lon)
head(lon)
lat <- ncvar_get(ncin, "lat", verbose = F)
nlat <- dim(lat)
head(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")
dim(tmp.array)
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")
nc_close(ncin)
# 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]
library(RColorBrewer)
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)
length(tmp.vec)
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)
length(tmp.vec.long)
tmp.mat <- matrix(tmp.vec.long, nrow = nlon * nlat, ncol = nt)
dim(tmp.mat)
head(na.omit(tmp.mat))
lonlat <- expand.grid(lon, lat)
tmp.df02 <- data.frame(cbind(lonlat, tmp.mat))
names(tmp.df02) <- c("lon","lat","pet_jan_2001",
"pet_feb_2001",
"pet_mar_2001",
"pet_apr_2001",
"pet_may_2001",
"pet_jun_2001",
"pet_jul_2001",
"pet_aug_2001",
"pet_sep_2001",
"pet_oct_2001",
"pet_nov_2001",
"pet_dec_2001",
"pet_jan_2002",
"pet_feb_2002",
"pet_mar_2002",
"pet_apr_2002",
"pet_may_2002",
"pet_jun_2002",
"pet_jul_2002",
"pet_aug_2002",
"pet_sep_2002",
"pet_oct_2002",
"pet_nov_2002",
"pet_dec_2002",
"pet_jan_2003",
"pet_feb_2003",
"pet_mar_2003",
"pet_apr_2003",
"pet_may_2003",
"pet_jun_2003",
"pet_jul_2003",
"pet_aug_2003",
"pet_sep_2003",
"pet_oct_2003",
"pet_nov_2003",
"pet_dec_2003",
"pet_jan_2004",
"pet_feb_2004",
"pet_mar_2004",
"pet_apr_2004",
"pet_may_2004",
"pet_jun_2004",
"pet_jul_2004",
"pet_aug_2004",
"pet_sep_2004",
"pet_oct_2004",
"pet_nov_2004",
"pet_dec_2004",
"pet_jan_2005",
"pet_feb_2005",
"pet_mar_2005",
"pet_apr_2005",
"pet_may_2005",
"pet_jun_2005",
"pet_jul_2005",
"pet_aug_2005",
"pet_sep_2005",
"pet_oct_2005",
"pet_nov_2005",
"pet_dec_2005",
"pet_jan_2006",
"pet_feb_2006",
"pet_mar_2006",
"pet_apr_2006",
"pet_may_2006",
"pet_jun_2006",
"pet_jul_2006",
"pet_aug_2006",
"pet_sep_2006",
"pet_oct_2006",
"pet_nov_2006",
"pet_dec_2006",
"pet_jan_2007",
"pet_feb_2007",
"pet_mar_2007",
"pet_apr_2007",
"pet_may_2007",
"pet_jun_2007",
"pet_jul_2007",
"pet_aug_2007",
"pet_sep_2007",
"pet_oct_2007",
"pet_nov_2007",
"pet_dec_2007",
"pet_jan_2008",
"pet_feb_2008",
"pet_mar_2008",
"pet_apr_2008",
"pet_may_2008",
"pet_jun_2008",
"pet_jul_2008",
"pet_aug_2008",
"pet_sep_2008",
"pet_oct_2008",
"pet_nov_2008",
"pet_dec_2008",
"pet_jan_2009",
"pet_feb_2009",
"pet_mar_2009",
"pet_apr_2009",
"pet_may_2009",
"pet_jun_2009",
"pet_jul_2009",
"pet_aug_2009",
"pet_sep_2009",
"pet_oct_2009",
"pet_nov_2009",
"pet_dec_2009",
"pet_jan_2010",
"pet_feb_2010",
"pet_mar_2010",
"pet_apr_2010",
"pet_may_2010",
"pet_jun_2010",
"pet_jul_2010",
"pet_aug_2010",
"pet_sep_2010",
"pet_oct_2010",
"pet_nov_2010",
"pet_dec_2010")
options(width = 110)
head(na.omit(tmp.df02, 20))
dim(na.omit(tmp.df02))
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
library(rgdal)
# 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")],
sep=""))
# 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 #
############################
library(rgdal)
nuts <- readOGR(".", layer = "NUTS_RG_60M_2006")
library(rgeos)
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
value.
# 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 https://stat.ethz.ch/mailman/listinfo/r-sig-geo We value your opinion. Please take a few minutes to share your comments
on
the service you received from the District by clicking on this link<
.
--
[image: Logo UHasselt] Mevrouw Janka Vanschoenwinkel
*Doctoraatsbursaal - PhD *
Milieueconomie - Environmental economics
T +32(0)11 26 86 96 | GSM +32(0)476 28 21 40
www.uhasselt.be/eec
Universiteit Hasselt | Campus Diepenbeek
Agoralaan Gebouw D | B-3590 Diepenbeek
Kantoor F11
Postadres: Universiteit Hasselt | Martelarenlaan 42 | B-3500 Hasselt
[image: Music For Life] Maak van UHasselt de #warmsteunief |
www.uhasselt.be/musicforlife
P Please consider the environment before printing this e-mail
[[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Dr. Michael Sumner Software and Database Engineer Australian Antarctic Division 203 Channel Highway Kingston Tasmania 7050 Australia [[alternative HTML version deleted]]