[R-sig-Geo] Translate a net-cdf file with different years to polygon regions.
Michael Sumner
mdsumner at gmail.com
Fri Apr 29 14:42:19 CEST 2016
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
> > <
> >
> https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_3.23/cruts.1506241137.v3.23/pet/
> > >
> > 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
> > <
> >
> https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_3.23/cruts.1506241137.v3.23/pet/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
> > <
> >
> http://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units/nuts
> > >
> > 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 :
> >
> >
> https://drive.google.com/folderview?id=0By9u5m3kxn9yfjZtdFZLcW82SWpzT1VwZXE1a3FtRGtSdEl1c1NvY205TGpack9xSFc2T2s&usp=sharing
> >
> > # 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<
> >
> http://my.sfwmd.gov/portal/page/portal/pg_grp_surveysystem/survey%20ext?pid=1653
> > >.
> >
>
>
>
> --
>
> [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]]
More information about the R-sig-Geo
mailing list