Re: [R-sig-Geo] Fwd: spatio temporal variogram
Hi! I was actually able to solve this using the spacetime and gstat packages in R. thanks for your help. regards, Jonesmus. On Wed, Feb 13, 2013 at 4:31 PM, Hodgess, Erin hodge...@uhd.edu wrote: Hi! When you make your shift from the first data set to the second, set it up such the the X1 - X4 are time series. For instance, if your data frame is called d1.df, d1.df$X1 - ts(d1.df$X1,start=1981,freq=12) (for monthly data) Hope this helps! Sincerely, Erin From: r-sig-geo-boun...@r-project.org [r-sig-geo-boun...@r-project.org] on behalf of Jonesmus Wambua [jonesmus.mu...@gmail.com] Sent: Wednesday, February 13, 2013 3:25 AM To: r-sig-geo@r-project.org Subject: [R-sig-Geo] Fwd: spatio temporal variogram Hi everyone, am trying to do temporal and cross correlation and also fit a spatio-temporal variogram with data in this format. | ID year y X_COORD Y_COORD 1 1 1981 60 36.98203 2.042057 2 2 1982 40 37.12720 2.016312 3 2 1983 20 37.12720 2.016312 4 3 1980 165 35.38416 1.714100 5 3 1981 26 35.38416 1.714100 6 3 1982 50 35.38416 1.714100| My challenge is getting the data in the right format. my years are only 5 but with over 500 cluster points. I have summarized the data by adding the data for each cluster per year so that my data now looks like below year X1 X2 X3 X4 1 1980 NA NA 3 NA 2 1981 1 NA 2 5 3 1982 NA 0 0 NA 4 1983 NA 0 NA NA 5 1984 NA NA NA NA where X1 are my clusters. when i tried to do a autocorrelation using the following code, acf(na.omit(som)) it works, but for cross-autocorrelation as follows acf(na.omit(som), xts) , it gives the following error Error in ts(x) : 'ts' object must have one or more observations can anyone help? thanks [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Need help with spatial data plotting
I suggest you run the code step by step and examine what happens. By reading help files along the way will bring you closer to understanding the code than you imagined. Here are a few hints I hope you may find helpful. data(eberg) This loads a dataset called `eberg` First, run coordinates(eberg). Error, right? Error in .checkNumericCoerce2double(obj) : cannot retrieve coordinates from non-numeric elements Now run coordinates(eberg) - ~X+Y This sets the default coordinates of this object. Inspect it with str(eberg) and you will see you're now dealing with a SpatialPointsDataFrame (very important to remember this!). When you imported the data, it was only a data.frame. It is now enhanced with coordinates, among other things. proj4string(eberg) - CRS(+init=epsg:31467) This will set the projection of those coordinates in a data.frame to espg:31467. eberg - eberg[runif(nrow(eberg)).2, ] First run nrow(eberg). This is the number of rows data.frame eberg has. Run runif(nrow(eberg)) and you will notice that you get a lot of values (nrow(eberg) many, as a matter of fact) between 0 and 1. By saying 2, you are asking, which values are smaller than 0.2? If a value is smaller than 0.2, it will be marked as TRUE and FALSE otherwise. By putting all this inside [..., ], it makes a convenient way of extracting random rows from a data.frame. Values before the comma means rows, after the comma means columns. Try eberg[1:5, 1:5]. See? If there are no commas, R assumes you mean columns. bubble(eberg[CLYMHT_A]) Calls a bubble plot (see ?bubble) and it takes column CLYMHT_A as data to draw bubbles. If this doesn't work for you, try the example provided in ?bubble: data(meuse) coordinates(meuse) - c(x, y) # promote to SpatialPointsDataFrame bubble(meuse, cadmium, maxsize = 2.5, main = cadmium concentrations (ppm), key.entries = 2^(-1:4)) plotKML(eberg[CLYMHT_A]) Ibidem, except you can now see the data displayed in Google Earth (or some other viewer). Cheers, Roman On Thu, Feb 14, 2013 at 10:06 AM, EPBaron ephr...@epbaron.com wrote: Thank you Roman. I had already looked at the plotKML tutorial. I found the examples fascinating, but I'm afraid it made me feel even more illiterate than when I started. In essence, it told me how to build an elegant watch when I'm just trying to tell what time it is. Or, to use the language analogy, it introduced me to poetry when all I can understand is simple sentences. The code that seemed most relevant to what I'm trying to do baffles me: data(eberg)coordinates(eberg) - ~X+Yproj4string(eberg) - CRS(+init=epsg:31467)eberg - eberg[runif(nrow(eberg)).2,]bubble(eberg[CLYMHT_A])plotKML(eberg[CLYMHT_A]) -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Need-help-with-spatial-data-plotting-tp7582619p7582626.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- In God we trust, all others bring data. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Spatio temporal variogram
Hi, i checked the mannuals and managed to sort out my problem. I just wanted to inquire about something that is not very clear in the specification of the variogram argument. what informs the specification of width and cutoff? am realizing that my results are changing depending on what i put there. som=variogram(PfPR~1, somal, width=10, cutoff=100) plot(som, ylab=time lag (years)) plot(som, map=FALSE, ylab=time lag (years)) Thanks in advance! [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Fwd: Spatio temporal variogram
Hi, i checked the mannuals and managed to sort out my problem. I just wanted to inquire about something that is not very clear in the specification of the variogram argument. what informs the specification of width and cutoff? am realizing that my results are changing depending on what i put there. i have 5 time points (in years) and 1164 spatial locations. som=variogram(PfPR~1, somal, width=10, cutoff=100) plot(som, ylab=time lag (years)) plot(som, map=FALSE, ylab=time lag (years)) -- Forwarded message -- From: Jonesmus Wambua jonesmus.mu...@gmail.com Date: Thu, Feb 14, 2013 at 2:08 PM Subject: Spatio temporal variogram To: r-sig-geo@r-project.org r-sig-geo@r-project.org Hi, i checked the mannuals and managed to sort out my problem. I just wanted to inquire about something that is not very clear in the specification of the variogram argument. what informs the specification of width and cutoff? am realizing that my results are changing depending on what i put there. som=variogram(PfPR~1, somal, width=10, cutoff=100) plot(som, ylab=time lag (years)) plot(som, map=FALSE, ylab=time lag (years)) Thanks in advance! [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Need help with spatial data plotting
Roman Luštrik wrote I suggest you run the code step by step and examine what happens. By reading help files along the way will bring you closer to understanding the code than you imagined. Here are a few hints I hope you may find helpful. ___ R-sig-Geo mailing list R-sig-Geo@ https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- In God we trust, all others bring data. ___ R-sig-Geo mailing list R-sig-Geo@ https://stat.ethz.ch/mailman/listinfo/r-sig-geo Roman, I know the rules of the SIG say I shouldn't simply reply to say 'thank you', but I really must. You're kindness is VERY much appreciated. I will work through the example again and carefully, and I'll try to replicate it with my data. br Wish me luck! -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Need-help-with-spatial-data-plotting-tp7582619p7582632.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Clip a contour with shapefile while using contourplot
Hi, I have done the interpolation for my data and I was able to create the contours in multipanel with the help of Pascal. Now, I want to clip the contour with the shapefile. I want only the portion of contour to be displayed which falls inside the boundary of the shapefile. The data mydata.csv can be found on https://www.dropbox.com/s/khi7nv0160hi68p/mydata.csv The data for shapefile can be found on https://www.dropbox.com/sh/ztvmibsslr9ocmc/YOtiwB8p9p THe code I have used so far is as follows: # Load Libraries library(latticeExtra) library(sp) library(rgdal) library(lattice) library(gridExtra) #Read Shapefile hello - readOGR(shape, layer=Export_Output_4) ## Project the shapefile to the UTM 16 zone proj4string(hello) - CRS(+proj=utm +zone=16 +ellps=WGS84) ## Read Contour data mydata - read.csv(mydata.csv) head(mydata ) summary(mydata) #Create a contourplot contourplot(Salinity ~ Eastings+Northings | Time, mydata, cuts=30,pretty=TRUE) Thank you so much. I would welcome any other ways to do this aside from contourplot and lattice. Best Regards, Janesh [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Change z coordinates of SpatialGridDataFrame
Dear Edzer, I transformed my SpatialGridDataFrame into a SpatiaPointsDataFrame (by the way is there any coercion method such as as.SpatialPointsDataFrame available in sp instead of creating a new one for I have not found it ?). But this does not change my problem: I need the vector of attributes to follow through the change in the z position. Right now the order of the values in the vector of attributes stays the same even if the z coordinates change. Here is my code starting from the summary of the initial SpatialGridDataFrame summary(tls.spgdf) Object of class SpatialGridDataFrame Coordinates: minmax [1,] 828099 828139 [2,] 101834 101874 [3,]322374 Is projected: TRUE proj4string : [+init=epsg:27561 +proj=lcc +lat_1=49.51 +lat_0=49.51 +lon_0=0 +k_0=0.999877341 +x_0=60 +y_0=20 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs] Grid attributes: cellcentre.offset cellsize cells.dim 1 828099.5140 2 101834.5140 3 322.5152 Data attributes: Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 0.00 0.00 0.007572 0.00 1.226000 # (i) 3D Matrix (x.coord - floor(coordinates(tls.spgdf)[,1])) (y.coord - floor(coordinates(tls.spgdf)[,2])) (z.coord - floor(coordinates(tls.spgdf)[,3])) tls.xyz2 - cbind(x.coord, y.coord, z.coord) # (ii) SpatialPoints tls.sp2 - SpatialPoints(tls.xyz2, proj4string = CRS(+init=epsg:27561)) # (iii) DataFrame (attribut z) tls.df2 - data.frame(d=tls.spgdf$d) # (iv) SpatialPointsDataFrame tls.spdf2 - SpatialPointsDataFrame(tls.xyz2, tls.df2) # Summary of DTM values summary(sge) Object of class SpatialGridDataFrame Coordinates: minmax s1 828099 828139 s2 101834 101874 Is projected: TRUE proj4string : [+init=epsg:27561 +proj=lcc +lat_1=49.51 +lat_0=49.51 +lon_0=0 +k_0=0.999877341 +x_0=60 +y_0=20 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs] Grid attributes: cellcentre.offset cellsize cells.dim s1 828099.5140 s2 101834.5140 Data attributes: Min. 1st Qu. MedianMean 3rd Qu.Max. 325.9 329.1 330.8 330.8 332.4 335.6 z.sge - round(sge[[1]]) z - coordinates(tls.spdf2)[,3] tls.spdf2$z - z nz - tls.spdf2[[2]] - z.sge[[1]] tls.spdf2$nz - nz summary(tls.spdf2) Object of class SpatialPointsDataFrame Coordinates: minmax x.coord 828099 828138 y.coord 101834 101873 z.coord322373 Is projected: NA proj4string : [NA] Number of points: 83200 Data attributes: d z nz Min. :0.00 Min. :322.0 Min. :-11.00 1st Qu.:0.00 1st Qu.:334.8 1st Qu.: 1.75 Median :0.00 Median :347.5 Median : 14.50 Mean :0.007572 Mean :347.5 Mean : 14.50 3rd Qu.:0.00 3rd Qu.:360.2 3rd Qu.: 27.25 Max. :1.226040 Max. :373.0 Max. : 40.00 (x.coord2 - coordinates(tls.spdf2)[,1]) (y.coord2 - coordinates(tls.spdf2)[,2]) (z.coord2 - tls.spdf2$nz) tls.xyz3 - cbind(x.coord2, y.coord2, z.coord2) tls.sp3 - SpatialPoints(tls.xyz3, proj4string = CRS(+init=epsg:27561)) tls.df2 - data.frame(d=tls.spdf2$d) tls.spdf3 - SpatialPointsDataFrame(tls.sp3, tls.df2) -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Change-z-coordinates-of-SpatialGridDataFrame-tp7582612p7582635.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] apply function to lists of rasters
This mapply problem has now been fixed for the next release of raster (= 2.1-6)
https://stat.ethz.ch/pipermail/r-devel/2013-February/065847.html
Also thanks to Roman Luštrik for off list suggestions. Robert
On Mon, Feb 11, 2013 at 9:14 AM, Robert J. Hijmans r.hijm...@gmail.com wrote:
mapply does indeed not seem to work here, further illustrated below. I do
not know why. However, in this case I think it is simpler not to use mapply
(see example at the very end).
library(raster)
Loading required package: sp
raster 2.1-3 (1-February-2013)
aa-matrix(c(2,1,3,4),nrow=2,ncol=2)
bb-matrix(c(5,7,4,6),nrow=2,ncol=2)
cc-matrix(c(12,13,17,16),nrow=2,ncol=2)
A-list(raster(aa),raster(t(aa)))
B-list(raster(bb),raster(t(bb)))
C-list(raster(cc),raster(t(cc)))
f1 - function(a,b,c) 1/(1-exp(-exp(-((a - b)/c
# this function does not work, I suppose the - function, as shortcut for
-1 * is not implemented for Raster objects
f1(A[[1]], B[[1]], C[[1]])
Error in -((a - b)/c) : invalid argument to unary operator
# You can make it work like this, because -1 * is understood
f2 - function(a,b,c) 1/(1-exp(-1*exp(-1*((a - b)/c
f2(A[[1]], B[[1]], C[[1]])
class : RasterLayer
dimensions : 2, 2, 4 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : in memory
names : layer
values : 1.257289, 1.529642 (min, max)
# now on to mapply
# simple fucnction
f3 - function(a,b,c) (a - b)/c
# works
D - mapply(f3, A,B,C)
# original function does not work
D - mapply(f2, A,B,C)
Error in as.character(sys.call(sys.parent())[[1]]) :
cannot coerce type 'closure' to vector of type 'character'
# probably because exp is part of the Math group generic and the
# method for Raster objects gets lost along the way.
# I need to investigate that. Ideas anyone?
# However, I think the good news is that a more rasteresque solution works:
f2(stack(A), stack(B), stack(C))
class : RasterBrick
dimensions : 2, 2, 4, 2 (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5 (x, y)
extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : in memory
names : layer.1, layer.2
min values : 1.257289, 1.257289
max values : 1.529642, 1.529642
# OR
overlay(stack(A), stack(B), stack(C), fun=f2)
class : RasterBrick
dimensions : 2, 2, 4, 2 (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5 (x, y)
extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : in memory
names : layer.1, layer.2
min values : 1.257289, 1.257289
max values : 1.529642, 1.529642
On Mon, Feb 11, 2013 at 5:31 AM, Roman Luštrik roman.lust...@gmail.com
wrote:
There may be a problem.
mapply(function(a, b, c) exp(a), A, B, C)
alas
exp(A[[1]])
class : RasterLayer
dimensions : 2, 2, 4 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : in memory
names : layer
values : 2.718282, 54.59815 (min, max)
Looks like exp() doesn't work from inside `mapply`. Robert?
# I changed c to d, because it's hell to debug `c`
D - mapply(function(a, b, d) {
browser()
1/(1 - exp(-exp(-((a - b)/d
}, A, B, C)
Browse[1] class(a)
[1] RasterLayer
attr(,package)
[1] raster
Browse[1] exp(a)
Error in as.character(sys.call(sys.parent())[[1]]) :
cannot coerce type 'closure' to vector of type 'character'
Cheers,
Roman
On Mon, Feb 11, 2013 at 2:14 PM, Lorenzo Alfieri
alfio...@hotmail.comwrote:
Dear list,
I'm trying to apply a function to lists of rasters, but without success.
I'm not sure if it feasible of if I should use a for loop.
Let's say a, b and c are three lists of argument to my function
1/(1-exp(-exp(-((a - b)/c (2 rasters each)
I'd like this function to be applied to these arguments, so that I get
as
result a list D of 2 rasters, where for example
D[[1]] = 1/(1-exp(-exp(-((A[[1]] - B[[1]])/C[[1]] and the same for
D[[2]]
I tryed with this but didn't work:
library(raster)
aa-matrix(c(2,1,3,4),nrow=2,ncol=2)
bb-matrix(c(5,7,4,6),nrow=2,ncol=2)
cc-matrix(c(12,13,17,16),nrow=2,ncol=2)
A-list(raster(aa),raster(t(aa)))
B-list(raster(bb),raster(t(bb)))
C-list(raster(cc),raster(t(cc)))
D-mapply(function(a,b,c) 1/(1-exp(-exp(-((a - b)/c, A,B,C)
Any suggestions?
Lorenzo
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
In God we trust, all others bring data.
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] measuring overlaps with many shapefiles
Hello, I am familiar with R and GIS separately, but not working with spatial objects in R. I have a PRIMARY shapefile and many SECONDARY shapefiles. All are polygons and have a string label column for features. I want to generate a nonspatial dataframe of overlap measurements with the following columns: 1) label of PRIMARY feature 2) SECONDARY shapefile name 3) label of SECONDARY feature 4) area of overlap in square miles I believe on way to proceed is to use v.to.rast and r.stats in spgrass6. But, I am not sure if that is the easiest way. Any suggestions would be appreciated. Thank you, -david ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] how to display a projected map on an rasterVis::layerplot?
Please advise me regarding overlay of an LCC-projected map on `raster` data plotted by `rasterVis::layerplot`. What I mean, why I ask: As discussed in detail @ http://stackoverflow.com/questions/14865507/how-to-display-a-projected-map-on-an-rlatticelayerplot (more detail than is feasible for an email post) I'm trying to overlay LCC-projected data from a netCDF file with an LCC-projected map. The netCDF uses the IOAPI convention (described in link), and also manages to hit yet again the `raster` filename-extension problem. (Dunno why, but that's twice in one week for me--bad karma?) The data's horizontal extents are CONUS, so I'm looking for an LCC CONUS map. In one example (in the link above), I use CRAN package=M3 to get (basically) a matrix of coordinates map.lines - M3::get.map.lines.M3.proj( file=lcc.file, database=state, units=m) and feed that to old-style graphics. But I'd much prefer to make this work like the other example, which resembles other code I have that is using `rasterVis::layerplot`. Unfortunately that other code is using unprojected/lon-lat data, and I must also handle projected (probably all LCC) data. TIA, Tom Roche tom_ro...@pobox.com ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
