Hello list I am stuck with a problem regarding the manipulation of two data sets. The first data set contains electromagnetic induction (EMI) readings of a field obtained with a high sampling density,14965 records in a field 1000 by 500 m. Each reading is geo-referenced. The second set contains soil related information (texture) taken from soil pits of the same field at a density of 75 x 75 m resulting in 111 readings. The position of each pit is geo-referenced as well. I would like to calculate an average of the EMI readings in the vicinity of the soil pits, say within a radius of 2 - 5 metres and relate this average to the soil texture data obtained from the pits. Averaging would also reduce the size of the EMI data set which at the moment is too big for ordinary kriging in GSTAT. Can this be done with a few lines of code? Regards _________________________________ Karl J Sommer, Department of Primary Industries, Catchment & Agriculture Services, PO Box 905 Mildura, VIC, Australia 3502 Tel: +61 (0)3 5051 4390 Fax +61 (0)3 5051 4534 Email: karl.sommer at dpi.vic.gov.au
data set manipulation, averaging
8 messages · karl.sommer at dpi.vic.gov.au, Edzer Pebesma, Roger Bivand +1 more
karl.sommer at dpi.vic.gov.au wrote:
Hello list I am stuck with a problem regarding the manipulation of two data sets. The first data set contains electromagnetic induction (EMI) readings of a field obtained with a high sampling density,14965 records in a field 1000 by 500 m. Each reading is geo-referenced. The second set contains soil related information (texture) taken from soil pits of the same field at a density of 75 x 75 m resulting in 111 readings. The position of each pit is geo-referenced as well. I would like to calculate an average of the EMI readings in the vicinity of the soil pits, say within a radius of 2 - 5 metres and relate this average to the soil texture data obtained from the pits. Averaging would also reduce the size of the EMI data set which at the moment is too big for ordinary kriging in GSTAT.
I assume you did not set a local neighbourhood for kriging? It seems to work for me, but takes indeed almost 2 Gb of RAM, > 15000^2*8 [1] 1.8e+09 bytes are needed for the single full covariance matrix alone.
Can this be done with a few lines of code?
use gstat with a maxdist (search radius) of 2 to 5 m, and a pure nugget variogram vgm(1, "Nug", 0) and you'll get the arithmetic mean of the data within the search radius. hth, -- Edzer
Regards
_________________________________ Karl J Sommer, Department of Primary Industries, Catchment & Agriculture Services, PO Box 905 Mildura, VIC, Australia 3502 Tel: +61 (0)3 5051 4390 Fax +61 (0)3 5051 4534 Email: karl.sommer at dpi.vic.gov.au _______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
7 days later
Hello, I'm sorry but I have lost track of how this is to be done. Any vector format would do, but I guess shapefiles are the most obvious. Is there a "spmaptools" package somewhere for doing this? I can't find it. Otherwise, how can I convert a "SpatialPolygons" object to a "polylist" (for writing by maptools)? Sorry if this is documented somewhere, but it's one of those things I had assumed could done and have only just attempted to find out how. Cheers, Mike.
On Wed, 30 Aug 2006, Michael Sumner wrote:
Hello, I'm sorry but I have lost track of how this is to be done. Any vector format would do, but I guess shapefiles are the most obvious. Is there a "spmaptools" package somewhere for doing this? I can't find it.
It is included in the maptools package now - the writePolyShape() function is the one you are looking for. If you also need to write a *.prj, use showWKT() in the rgdal package. Code is getting ahead of documentation here ... Roger
Otherwise, how can I convert a "SpatialPolygons" object to a "polylist" (for writing by maptools)? Sorry if this is documented somewhere, but it's one of those things I had assumed could done and have only just attempted to find out how. Cheers, Mike.
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no
It is included in the maptools package now - the writePolyShape() function is the one you are looking for. If you also need to write a *.prj, use
Ah, thanks. But how do I convert an "SpatialPolygons" object, or even the @polygons list to a "polylist"? Sorry . . . an example? Cheers, Mike.
On Thu, 31 Aug 2006, Michael Sumner wrote:
It is included in the maptools package now - the writePolyShape() function is the one you are looking for. If you also need to write a *.prj, use
Ah, thanks. But how do I convert an "SpatialPolygons" object, or even the @polygons list to a "polylist"? Sorry . . . an example?
What is a shapefile? It has a *.shp, a *.shx, _and_ a *.dbf at least. So
you can't go from a SpatialPolygons object to a shapefile, it has to be
SpatialPolygonsDataFrame, even if the only column in the data frame is
just a holder. So create a data frame (with row names equal to the
Polygons ID values), make a SpatialPolygonsDataFrame, and
writePolyShape().
library(maptools)
xx <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
SPxx <- as(xx, "SpatialPolygons")
class(SPxx)
IDs <- sapply(slot(SPxx, "polygons"), function(x) slot(x, "ID"))
df <- data.frame(rep(0, length(IDs)), row.names=IDs)
SPDFxx <- SpatialPolygonsDataFrame(SPxx, df)
class(SPDFxx)
summary(SPDFxx)
tf <- tempfile()
writePolyShape(SPDFxx, tf)
getinfo.shape(tf)
The conversion stuff is probably irrelevant, it is internal and subject to
change. It's in the source if you do need it.
Roger
Cheers, Mike.
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no
What is a shapefile? It has a *.shp, a *.shx, _and_ a *.dbf at least. So you can't go from a SpatialPolygons object to a shapefile, it has to be SpatialPolygonsDataFrame, even if the only column in the data frame is just a holder. So create a data frame (with row names equal to the Polygons ID values), make a SpatialPolygonsDataFrame, and writePolyShape().
Oh sorry, I had missed that "writePolyShape" function. It's aliased to its read counterpart, and I thought I was seeing write.polylistShape. Thanks!
On Thu, 31 Aug 2006, Michael Sumner wrote:
What is a shapefile? It has a *.shp, a *.shx, _and_ a *.dbf at least. So you can't go from a SpatialPolygons object to a shapefile, it has to be SpatialPolygonsDataFrame, even if the only column in the data frame is just a holder. So create a data frame (with row names equal to the Polygons ID values), make a SpatialPolygonsDataFrame, and writePolyShape().
Oh sorry, I had missed that "writePolyShape" function. It's aliased to its read counterpart, and I thought I was seeing write.polylistShape.
I'll split it out if you like - in the help.start() browser based help system, each alias is listed separately, so you would have seen it there.
Thanks!
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no