Thanks for explaining this in such detail! I have a greater appreciation for the importance of thinking about these topological issues and for the role machine precision plays.
I constructed a series of simple examples to demonstrate to myself how these sorts of problems arise and how setScale, set_RGEOS_polyThreshold, set_RGEOS_dropSlivers, and set_RGEOS_warnSlivers work. In case someone else stumbles on this thread looking for similar clarification, I've pasted these examples below, and they also appear in this gist: https://gist.github.com/mstrimas/1b4a4b93a9d4a158bce4 library(sp) library(raster) library(rgeos) set_RGEOS_polyThreshold(0) set_RGEOS_warnSlivers(TRUE) set_RGEOS_dropSlivers(FALSE) # function to rigidly shift polygon # only works for simple polygon objects with a single polygon shift_poly <- function(p, x, y) { p@polygons[[1]]@Polygons[[1]]@coords[,'x'] <- p@polygons[[1]]@Polygons[[1]]@coords[,'x'] + x p@polygons[[1]]@Polygons[[1]]@coords[,'y'] <- p@polygons[[1]]@Polygons[[1]]@coords[,'y'] + y row.names(p) <- paste0('shifted', row.names(p)) return(p) } p1 <- readWKT("POLYGON((0 0,0 1,1 1,1 0,0 0))") row.names(p1) <- '1' p2 <- readWKT("POLYGON((1 0,1 1,2 1,2 0,1 0))") row.names(p2) <- '2' plot(rbind(p1, p2)) plot(shift_poly(p2, -0.5, 0.1), border='orange', add=T, lty=2, lwd=1) # default scale of 10^8 # behaviour of topology operations depends on scale! setScale(1e8) # shift horizontally by small amount so no longer touching plot(gUnion(p1, shift_poly(p2, 1e-1, 0))) plot(gUnion(p1, shift_poly(p2, 1e-8, 0))) plot(gUnion(p1, shift_poly(p2, 1e-9, 0))) gIntersects(p1, p2) gIntersects(p1, shift_poly(p2, 1e-1, 0)) gIntersects(p1, shift_poly(p2, 1e-8, 0)) gIntersects(p1, shift_poly(p2, 1e-9, 0)) # polygons with coordinates different by less that the scale set by setScale() # are considered to intersect and are merge together by union operations # p1 and p2 share a boundary exactly so the intersection of the 2 is a line class(gIntersection(p1, p2)) # shift right polygon horizontally slightly to it is just overlapping or just # separated from the left polygon and consider the results gIntersection(p1, shift_poly(p2, -1e-9, 0)) # overlap > scale => polygon gIntersection(p1, shift_poly(p2, -1e-8, 0)) # overlap < scale => line gIntersection(p1, p2) # exactly touching => line gIntersection(p1, shift_poly(p2, 1e-9, 0)) # separation < scale => line gIntersection(p1, shift_poly(p2, 1e-8, 0)) # separation > scale => NULL # lower scale to 10^4 setScale(1e4) plot(gUnion(p1, shift_poly(p2, 1e-1, 0))) plot(gUnion(p1, shift_poly(p2, 1e-4, 0))) plot(gUnion(p1, shift_poly(p2, 1e-5, 0))) gIntersects(p1, p2) gIntersects(p1, shift_poly(p2, 1e-1, 0)) gIntersects(p1, shift_poly(p2, 1e-4, 0)) gIntersects(p1, shift_poly(p2, 1e-5, 0)) gIntersection(p1, shift_poly(p2, -1e-4, 0)) # overlap > scale => polygon gIntersection(p1, shift_poly(p2, -1e-5, 0)) # overlap < scale => line gIntersection(p1, p2) # exactly touching => line gIntersection(p1, shift_poly(p2, 1e-5, 0)) # separation < scale => line gIntersection(p1, shift_poly(p2, 1e-4, 0)) # separation > scale => NULL # consider the effect of setting polyThreshold and turning on sliver warning # this will identify slivers resulting from topology operations setScale(1e4) set_RGEOS_polyThreshold(1e-2) set_RGEOS_warnSlivers(TRUE) # shift horizontally by increasing amount so just overlapping gi <- gIntersection(p1, shift_poly(p2, -1e-5, 0)) class(gi) # overlap < scale => line, i.e. assumes just touching # warnings raised in next 2 cases because: # a. overlap > scale => polygon on intersection # b. resulting polygon area < polyThreshold => sliver gIntersection(p1, shift_poly(p2, -1e-4, 0)) gIntersection(p1, shift_poly(p2, -1e-3, 0)) # warnings raised in next 2 cases because: # a. overlap > scale => polygon on intersection # b. resulting polygon area > polyThreshold => not a sliver gIntersection(p1, shift_poly(p2, -1e-2, 0)) gIntersection(p1, shift_poly(p2, -1e-1, 0)) # with a lower threshold set_RGEOS_polyThreshold(1e-3) # this still raises a warning gIntersection(p1, shift_poly(p2, -1e-4, 0)) # but this doesn't since resulting polygon is at the threshold now gIntersection(p1, shift_poly(p2, -1e-3, 0)) # not that it isn't the linear overlap that triggers the warning, it is that the # area of the resulting polygons are below the threshold # no warning gi1 <- gIntersection(p1, shift_poly(p2, -1e-3, 0)) # now a warning is raised because a slight shift in the y direction has caused # the polygons the resulting polygon to be just less than the 1e-3 threshold gi2 <- gIntersection(p1, shift_poly(p2, -1e-3, 1e-3)) gArea(gi1) / get_RGEOS_polyThreshold() gArea(gi2) / get_RGEOS_polyThreshold() # rgeos can also be set to automatically remove these slivers class(gIntersection(p1, shift_poly(p2, -1e-4, 0))) # SpaitalPolygons set_RGEOS_dropSlivers(TRUE) class(gIntersection(p1, shift_poly(p2, -1e-4, 0))) # NULL set_RGEOS_dropSlivers(FALSE) # slivers can also be generated as a result of union operations setScale(1e4) set_RGEOS_polyThreshold(1e-2) set_RGEOS_warnSlivers(TRUE) set_RGEOS_dropSlivers(FALSE) p3 <- readWKT("POLYGON((0 1,0 2,2 2,2 1,0 1))") row.names(p3) <- '3' p4 <- readWKT("POLYGON((0 -1,0 0,2 0,2 -1,0 -1))") row.names(p4) <- '4' plot(rbind(p1, p3, p4)) plot(p2, add=T, col='red') # offset the middle right (i.e. red) square slightly to the right # not picked up since middle edge misalignment is < scale pp <- rbind(p1, p3, p4, shift_poly(p2, 1e-5, 0)) plot(gUnaryUnion(pp)) # misalignment picked up since = scale => a very narrow hole in center # warning raised because hole area is < polyThreshold pp <- rbind(p1, p3, p4, shift_poly(p2, 1e-4, 0)) plot(gUnaryUnion(pp)) # misalignment picked up since = scale => a very narrow hole in center # warning not raised because hole area is > polyThreshold pp <- rbind(p1, p3, p4, shift_poly(p2, 1e-2, 0)) plot(gUnaryUnion(pp)) # the fact that this is a hole and not a vertical line becomes apparent when # the shift is bigger pp <- rbind(p1, p3, p4, shift_poly(p2, 1e-1, 0)) plot(gUnaryUnion(pp)) # finally, set_RGEOS_dropSlivers can be used to remove these slivers # and fix the topology set_RGEOS_dropSlivers(TRUE) pp <- rbind(p1, p3, p4, shift_poly(p2, 1e-4, 0)) plot(gUnaryUnion(pp)) set_RGEOS_dropSlivers(FALSE) On Thu, Nov 5, 2015 at 2:49 AM, Roger Bivand <roger.biv...@nhh.no> wrote: > On Thu, 5 Nov 2015, Roger Bivand wrote: > >> On Wed, 4 Nov 2015, Matt Strimas-Mackey wrote: >> >>> The original study area shapefile is a boundary of the Indonesia half >>> of New Guinea. The file as well as the code to construct the hexagonal >>> grids are here: >>> >>> https://www.dropbox.com/sh/ff8v08p3ambqcbs/AAAPBlGP4fthdmZhrto7oIuCa?dl=0 >>> >>> Since it's a large area, generating the grid takes a long time, so >>> I've also included code for a small subset of the original >>> shapefile--one small offshore island. >> >> >> An equivalent scaling fix is to change the coordinate units from metres to >> kilometres: >> >> library(rgdal) >> study_area2 <- spTransform(study_area, CRS("+proj=aea +lat_1=-7.42 >> +lat_2=-0.62 +lat_0=-9.119 +lon_0=128.91 +x_0=0 y_0=0 +datum=WGS84 >> +units=km +no_defs +ellps=WGS84 +towgs84=0,0,0")) >> size <- sqrt(2 * 1e2 / sqrt(3)) # reduce to suit >> set.seed(1) >> hex_points <- spsample(study_area2[2,], type="hexagonal", cellsize=size) >> hex_grid <- HexPoints2SpatialPolygons(hex_points) >> hex_union <- gUnaryUnion(hex_grid) >> plot(hex_grid, col='lightgrey') >> plot(hex_union, border='orange', lwd=3, add=T) >> >> An advantage of these scaling issues being made visible is that we see >> that computers really do not operate in continuous space, and that >> computational geometry actually matters, I suppose. >> >> I'll check the underlying code generating the hexagons to see why the >> nodes (where hexagon boundaries meet) appear to slide apart at the edge of >> machine precision. > > > A potentially robust approach at default scale uses the dx= argument to > HexPoints2SpatialPolygons(): > > hex_grid <- HexPoints2SpatialPolygons(hex_points, dx=size) > > Setting: > > set_RGEOS_polyThreshold(1e-2) # for example > set_RGEOS_warnSlivers(TRUE) > > shows the remaining slivers, and: > > set_RGEOS_dropSlivers(TRUE) > > drops them. It is still possible that an inward dangle will get through > (zero area line in from boundary point, but part of the single boundary). > > Setting dx= to the same value as cellsize= in the spsample call seems to be > crucial, avoiding an approximation marked in the code as: > > # EJP; changed: > # how to figure out dx from a grid? THK suggested: > #dx <- hexGrid$x[2] - hexGrid$x[1] > # and the following will also not always work: > > in sp:::genPolyList(), so there was a warning there against taking the > default. > > Roger > > >> >> Roger >> >> >>> >>> Finally, some more odd behaviour. I noticed each time I run this code, >>> the dissolve mistakes change, i.e. different boundaries are >>> erroneously kept. However, using set.seed() makes the errors the same >>> each time for a given seed, and changing the seed yields a different >>> set of errors. Example in the code in the dropbox link and copied >>> here: >>> >>> library(sp) >>> library(raster) >>> library(rgeos) >>> >>> # just a subset of full shapefile >>> set.seed(1) >>> size <- sqrt(2 * 1e8 / sqrt(3)) >>> study_area <- shapefile('papua.shp') >>> hex_points <- spsample(study_area[2,], type="hexagonal", cellsize=size) >>> hex_grid <- HexPoints2SpatialPolygons(hex_points) >>> hex_union <- gUnaryUnion(hex_grid) >>> plot(hex_grid, col='lightgrey') >>> plot(hex_union, border='orange', lwd=3, add=T) >>> >>> # results chage according to seed >>> set.seed(100) >>> size <- sqrt(2 * 1e8 / sqrt(3)) >>> study_area <- shapefile('papua.shp') >>> hex_points <- spsample(study_area[2,], type="hexagonal", cellsize=size) >>> hex_grid <- HexPoints2SpatialPolygons(hex_points) >>> hex_union <- gUnaryUnion(hex_grid) >>> plot(hex_grid, col='lightgrey') >>> plot(hex_union, border='orange', lwd=3, add=T) >>> >>> On Wed, Nov 4, 2015 at 8:57 AM, Roger Bivand <roger.biv...@nhh.no> >>> wrote: >>> > On Wed, 4 Nov 2015, Matt Strimas-Mackey wrote: >>> > > > Thanks, lots of useful info here. I've never seen the setScale() >>> > > function; I don't think it's mentioned in the gUnaryUnion help. This >>> > > saves me a lot of headache! >>> > > > > For what it's worth, the invalid geometry is an artifact of the >>> > > reproducible example I created. The original hexagonal grid is >>> > > produced with >>> > > g <- spsample(study_area, type="hexagonal", cellsize=size) >>> > > hex_grid <- HexPoints2SpatialPolygons(g) >>> > > > OK, thanks - this is useful. Could you please make available the >>> > > > study > area >>> > object in some way - so that we can re-create g and see how >>> > HexPoints2SpatialPolygons() creates the artefact Mike spotted >>> > (although > this >>> > looks like numeric fuzz - 'changing "1287248.96712942" to >>> > "1287248.96712943"' is a change in the 15th digit, which is on the > >>> > precision >>> > edge of the "double" storage mode. If we can revisit functions >>> > creating >>> > SpatialPolygons objects to ensure that they are GEOS-compatible for >>> > the >>> > default scale of 1e+8, we'll be more secure. >>> > > Roger >>> > > > > > > And this object passes gIsValid() and clgeo_GeometryReport() >>> > > > > > > without >>> > > any problems, yet still has the dissolving issue. Regardless, all is >>> > > solved with setScale(). >>> > > > > Thanks! >>> > > > > M >>> > > > > On Wed, Nov 4, 2015 at 6:30 AM, Michael Sumner >>> > > > > <mdsum...@gmail.com> > > wrote: >>> > > > > > > > > > > > > On Thu, 5 Nov 2015 at 01:10 Roger Bivand >>> > > > > > > > > > > > > <roger.biv...@nhh.no> > > > wrote: >>> > > > > > > > > > > > > On Wed, 4 Nov 2015, Michael Sumner wrote: >>> > > > > > > > > > Thanks for all this detail Roger, is there a way to >>> > > > > > > > > > "re-build" a >>> > > > > > spatial >>> > > > > > object so that the given scale setting is applied? Are there >>> > > > > > any >>> > > > > > general >>> > > > > > rounding or "orthogonalize" functions in the Spatial suite? >>> > > > > > > > > > > > > No, not really. In this case, the very detailed >>> > > > > > > > > > > > > coordinate > > > > measurements >>> > > > > may have made things worse, or possibly using Polygon rather >>> > > > > than >>> > > > > Polygons >>> > > > > objects, or not building the Polygons object with a proper >>> > > > > comment >>> > > > > attribute, I don't know. rgeos::gIsValid(spoly) is FALSE, and >>> > > > > > > > > > > > > >>> > > > > > > > > > > > > cleangeo::clgeo_GeometryReport(clgeo_CollectionReport(clgeo_Clean(spoly))) >>> > > > > > > > > says the same (based on rgeos). I suggest working with >>> > > > > > > > > Emmanuel > > > > Blondel >>> > > > > (cleangeo maintainer) to extend cleangeo. GEOS is looking at > > >>> > > > > > > allowing >>> > > > > users to manipulate precision models, not just scale, but I'm > >>> > > > > > > > uncertain >>> > > > > about that. >>> > > > > > > > > Running spoly into GRASS and back out (GRASS builds >>> > > > > > > > > topology on > > > > import) >>> > > > > shows a different error, the object seems to be problematic. >>> > > > > > > > > > > It can be fixed by changing "1287248.96712942" to > >>> > > > > > > > > > > > > "1287248.96712943", so >>> > > > really the creator should not have had those values on input. > > >>> > > > > There's no >>> > > > easy way out without topology. >>> > > > > > > This brought out some interesting issues for work I've been >>> > > > > > > doing > > > using >>> > > > the >>> > > > "near-Delaunay triangulation" in RTriangle, and that requires a >>> > > > normalized >>> > > > set of vertices (no duplicate vertices) which on its own presents >>> > > > interesting problems. I have a related issue where a parallel > > >>> > > > > (latitude) >>> > > > line needs many vertices to look "smooth" on a polar projection, >>> > > > but > > > when >>> > > > building a mesh with triangles it's really best to allow >>> > > > relatively >>> > > > coarse >>> > > > segmented boundaries rather than have many elements at parallels. >>> > > > > > > The >>> > > > Triangle library does not consider these hexagon coordinates to be >>> > > > duplicates, so there are two vertical segments between the two > > >>> > > > > bottom >>> > > > polys >>> > > > at >>> > > > > > > points(coordinates(as(as(spoly, "SpatialLines"), > > > >>> > > > > > > "SpatialPoints"))[c(3, >>> > > > 4), >>> > > > ] ) > > > > > > Thanks for the reminder about cleangeo, I'll have >>> > > > closer look. >>> > > > > > > cheers, Mike. >>> > > > > > > > > > > > Roger >>> > > > > > > > > > > > > > > Cheers, Mike. >>> > > > > > > > > > > On Wed, 4 Nov 2015 at 18:16 Roger Bivand >>> > > > > > > > > > > <roger.biv...@nhh.no> > > > > > wrote: >>> > > > > > > > > > > > On Tue, 3 Nov 2015, Matt Strimas-Mackey wrote: >>> > > > > > > > > > > > > > I'm working with a regular hexagonal grid >>> > > > > > > > > > > > > > stored as SPDF. At > > > > > > > some >>> > > > > > > > point I subset this SPDF, then want to combine all >>> > > > > > > > adjacent > > > > > > > hexagons >>> > > > > > > > together so that each contiguous set of hexagons is a >>> > > > > > > > single > > > > > > > polygon. >>> > > > > > > > I'm doing this last step using gUnaryUnion (or > > > > > > >>> > > > > > > > > gUnionCascaded, not >>> > > > > > > > clear what the different is actually). The problem is that >>> > > > > > > > > > > > > > > in some >>> > > > > > > > cases boundaries between clearly adjacent polygons are not >>> > > > > > > > > > > > > > > dissolved. >>> > > > > > > > > > > > > > > Example: >>> > > > > > > > > > > > > > > ## Create three adjacent hexagons >>> > > > > > > > library(sp) >>> > > > > > > > library(rgeos) >>> > > > > > > > p1 <- Polygon(cbind( >>> > > > > > > > c(1276503.26781119, 1281876.11747031, >>> > > > > > > > 1287248.96712942, >>> > > > > > > > 1287248.96712942, 1281876.11747031, >>> > > > > > > > 1276503.26781119, >>> > > > > > > > > > > > > 1276503.26781119), >>> > > > > > > > > > > > > > > c(204391.40834643, 207493.42454344, >>> > > > > > > > > > > > > > > 204391.40834643, >>> > > > > > > > > > > > > 198187.37595242, >>> > > > > > > > > > > > > > > 195085.35975541, 198187.37595242, >>> > > > > > > > > > > > > > > 204391.40834643))) >>> > > > > > > > p2 <- Polygon(cbind( >>> > > > > > > > c(1287248.96712943, 1292621.81678854, >>> > > > > > > > 1297994.66644766, >>> > > > > > > > 1297994.66644766, 1292621.81678854, >>> > > > > > > > 1287248.96712943, >>> > > > > > > > > > > > > 1287248.96712943), >>> > > > > > > > > > > > > > > c(204391.40834643, 207493.42454344, >>> > > > > > > > > > > > > > > 204391.40834643, >>> > > > > > > > > > > > > 198187.37595242, >>> > > > > > > > > > > > > > > 195085.35975541, 198187.37595242, >>> > > > > > > > > > > > > > > 204391.40834643))) >>> > > > > > > > p3 <- Polygon(cbind( >>> > > > > > > > c(1281876.11747031, 1287248.96712943, >>> > > > > > > > 1292621.81678854, >>> > > > > > > > 1292621.81678854, 1287248.96712943, >>> > > > > > > > 1281876.11747031, >>> > > > > > > > > > > > > 1281876.11747031), >>> > > > > > > > > > > > > > > c(213697.45693745, 216799.47313446, >>> > > > > > > > > > > > > > > 213697.45693745, >>> > > > > > > > > > > > > 207493.42454344, >>> > > > > > > > > > > > > > > 204391.40834643, 207493.42454344, >>> > > > > > > > > > > > > > > 213697.45693745))) >>> > > > > > > > spoly <- SpatialPolygons(list(Polygons(list(p1, p2, p3), > >>> > > > > > > > > > > > > > 's1'))) >>> > > > > > > > plot(gUnaryUnion(spoly)) >>> > > > > > > > > > > > > > > > > > > No, this is just numerical fuzz: >>> > > > > > > > > > > > > plot(spoly) >>> > > > > > > plot(gUnaryUnion(spoly), border="green", lty=3, lwd=2, > > > >>> > > > > > > > > > add=TRUE) >>> > > > > > > oS <- getScale() >>> > > > > > > # default 1e+8 >>> > > > > > > setScale(1e+4) >>> > > > > > > plot(gUnaryUnion(spoly), border="orange", lwd=2, add=TRUE) >>> > > > > > > setScale(oS) >>> > > > > > > > > > > > > JTS, GEOS, and consequently rgeos by default >>> > > > > > > > > > > > > shift all > > > > > > coordinates to >>> > > > > > > an >>> > > > > > > integer grid after multiplying by a scale factor (finding > >>> > > > > > > > > > > > integer >>> > > > > > > matches >>> > > > > > > is much easier than real matches). If the scaling is too > > >>> > > > > > > > > > > detailed (in >>> > > > > > > some >>> > > > > > > cases), the operations do not give the expected outcomes. >>> > > > > > > > > > > > > There is work in progress in GEOS and JTS to >>> > > > > > > > > > > > > provide other > > > > > > scaling >>> > > > > > > options >>> > > > > > > and models, and to permit iteration over scaling values >>> > > > > > > until > > > > > > a >>> > > > > > > "clean" >>> > > > > > > result is obtained (for some meanings of clean). >>> > > > > > > > > > > > > gUnionCascaded() was the only possible function >>> > > > > > > > > > > > > for GEOS < > > > > > > 3.3.0, from >>> > > > > > > GEOS 3.3.0 gUnaryUnion() is available and the prefered and > >>> > > > > > > > > > > > more >>> > > > > > > efficient >>> > > > > > > route. This is explained on the help page. >>> > > > > > > > > > > > > Hope this clarifies, >>> > > > > > > > > > > > > Roger >>> > > > > > > > > > > > > > > > > > > > > Note that p2 and p3 are >>> > > > > > > > > > > > > > > > > > > > > dissolved together, but p1 is > > > >>> > > > > > > > > > > > > > > > > > > > > > > > > separate. The >>> > > > > > > > shared edge of p1 and p2 is: >>> > > > > > > > p1: >>> > > > > > > > [2,] 1281876 207493.4 >>> > > > > > > > [3,] 1287249 204391.4 >>> > > > > > > > p2: >>> > > > > > > > [5,] 1287249 204391.4 >>> > > > > > > > [6,] 1281876 207493.4 >>> > > > > > > > > > > > > > > So, exactly the same apart from the order. I >>> > > > > > > > > > > > > > > originally > > > > > > > thought this >>> > > > > > > > difference in order might be the problem, but this doesn't >>> > > > > > > > > > > > > > > seem to be >>> > > > > > > > an issue with in this example, where order is also >>> > > > > > > > flipped: >>> > > > > > > > sss <- rasterToPolygons(raster(nrow=2, ncol=2, xmn=0, >>> > > > > > > > xmx=2, > > > > > > > ymn=0, >>> > > > > > > > ymx=2, vals=1:4)) >>> > > > > > > > lapply(sss@polygons, function(x) slot(x, > > > > > > > >>> > > > > > > > 'Polygons')[[1]]@coords) >>> > > > > > > > plot(sss) >>> > > > > > > > plot(gUnaryUnion(sss)) >>> > > > > > > > > > > > > > > Session Info: >>> > > > > > > > R version 3.2.2 (2015-08-14) >>> > > > > > > > Platform: x86_64-apple-darwin13.4.0 (64-bit) >>> > > > > > > > Running under: OS X 10.10.5 (Yosemite) >>> > > > > > > > > > > > > > > Message when rgeos is loaded: >>> > > > > > > > > > > > > > > rgeos version: 0.3-14, (SVN revision 511) >>> > > > > > > > GEOS runtime version: 3.5.0-CAPI-1.9.0 r0 >>> > > > > > > > Linking to sp version: 1.2-0 >>> > > > > > > > Polygon checking: TRUE >>> > > > > > > > > > > > > > > Any help on how to get these polygons to >>> > > > > > > > > > > > > > > dissolve is > > > > > > > appreciated. >>> > > > > > > > > > > > > > > M >>> > > > > > > > > > > > > > > >>> > > > > > > > > > > > > > > _______________________________________________ >>> > > > > > > > R-sig-Geo mailing list >>> > > > > > > > R-sig-Geo@r-project.org >>> > > > > > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> > > > > > > > > > > > > > > > > > > > -- >>> > > > > > > Roger Bivand >>> > > > > > > Department of Economics, Norwegian School of Economics, >>> > > > > > > Helleveien 30, N-5045 Bergen, Norway. >>> > > > > > > voice: +47 55 95 93 55; fax +47 55 95 91 00 >>> > > > > > > e-mail: roger.biv...@nhh.no >>> > > > > > > > > > > > > _______________________________________________ >>> > > > > > > R-sig-Geo mailing list >>> > > > > > > R-sig-Geo@r-project.org >>> > > > > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> > > > > > > > > > > > > > > > > > > > -- >>> > > > > Roger Bivand >>> > > > > Department of Economics, Norwegian School of Economics, >>> > > > > Helleveien 30, N-5045 Bergen, Norway. >>> > > > > voice: +47 55 95 93 55; fax +47 55 95 91 00 >>> > > > > e-mail: roger.biv...@nhh.no >>> > > > > > > > > > > > -- >>> > Roger Bivand >>> > Department of Economics, Norwegian School of Economics, >>> > Helleveien 30, N-5045 Bergen, Norway. >>> > voice: +47 55 95 93 55; fax +47 55 95 91 00 >>> > e-mail: roger.biv...@nhh.no >>> > >> >> >> > > -- > Roger Bivand > Department of Economics, Norwegian School of Economics, > Helleveien 30, N-5045 Bergen, Norway. > voice: +47 55 95 93 55; fax +47 55 95 91 00 > e-mail: roger.biv...@nhh.no > _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo