Re: [R-sig-Geo] Cleaning small spatial polygons

[Roger Bivand]

On Wed, 4 Nov 2015, Roger Bivand wrote:


On Wed, 4 Nov 2015, Eduardo Diez wrote:


 Is there any way of doing this or should i forget it and go on using GRASS
 through rgrass7?


I suggest working with Emmanuel Blondel (cleangeo maintainer) to extend 
cleangeo (also suggested in an earlier thread today, with apologies to 
Emmanuel for picking on him!). That package already uses rgeos, and is a 
logical place to put GRASS v.clean-like functionality (maybe even 
encapsulating using GRASs via rgrass7 and a throwaway location).


At the moment it is messy, though some rgeos internals do do something like 
this, but it isn't exposed to users.


A possibility is to use:

set_RGEOS_polyThreshold(1e-2) # for example
set_RGEOS_warnSlivers(TRUE)

shows the remaining slivers, and:

set_RGEOS_dropSlivers(TRUE)

drops them when using for example:


t1 <- gBuffer(, byid=TRUE, width=0)


A buffer of zero width should not have side-effects, but your mileage may 
vary. It will only remove slivers, not dangles. I haven't tried it when it 
also finds a Polygons object under the threshold, which was your initial 
problem.


Roger



Roger



 Thanks

 2015-10-19 15:03 GMT-03:00 Eduardo Diez :

>  Ok. So here's a link to a zip file that contains two shapefiles:
>  - pol_to_be_cleaned: the layer from which i'd like to remove small
>  polygons
>   - pol_cleaned: the layer cleaned with the function v.clean rmarea
> 
>  http://1drv.ms/1GmRWS7
> 
>  The threshold i used for cleaning was 3000 (meaning 3000 squared 
>  meters).
> 
>  Although i do project it before sending it to GRASS, according to the

>  official help page it should be able to handle it:
>  "Threshold must always be in square meters, also for latitude-longitude
>  locations or locations with units other than meters"
> 
>  Thanks
> 
>  2015-10-19 8:28 GMT-03:00 Roger Bivand :
> 
> >  On Fri, 16 Oct 2015, Eduardo Diez wrote:
> > 
> >  Dear list,
> > >  I'm willing to know if any knows a way of performing tha same thing 
> > >  i'm
> > >  doing through rgrass7 with GRASS when I execute the function v.clean 
> > >  with

> > >  "rmarea" as the tool argument. That is:
> > > 
> > >  "The rmarea tool removes all areas <= thresh. The longest boundary 
> > >  with

> > >  an
> > >  adjacent area is removed or all boundaries if there is no adjacent 
> > >  area.

> > >  Area categories are not combined when a small area is merged with a
> > >  larger
> > >  area."
> > > 
> > >  Basically i have raster of zones within a field. I convert it to

> > >  SpatialPolygonsDataFrame and in order to leave only the more
> > >  important/meaningful ones i remove the small/sliver with this tool. 
> > >  In

> > >  general it works fine but having to call an external software with a
> > >  specific version makes the script less portable and you have to be
> > >  careful
> > >  with updates and such. Also you have to write rasters and shapefiles 
> > >  back

> > >  and forth as GRASS can't work with in-memory objects.
> > > 
> > 
> >  Could you please provide an example of a built-in or contributed data 
> >  set
> >  (URL, not attachment) with the slivers you mention, so that we know 
> >  that we

> >  are addressing your problem? I don't think that:
> > 
> >  https://cran.r-project.org/web/packages/cleangeo/index.html
> > 
> >  does this, as it seems to try to repair broken geometries.
> > 
> >  Also note that you need to specify that the area threshold is in a 
> >  square

> >  planar metric - dropping slivers in unprojected geometries may be more
> >  complicated.
> > 
> >  Roger
> > 
> > 
> > >  Does someone know a way of doing this in plain R?
> > > 
> > >  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
> > > 
> > > 
> >  --

> >  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

Re: [R-sig-Geo] gUnaryUnion Not Dissolving Correctly

[Roger Bivand]

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  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  
> >  wrote:
> > > 
> > > 
> > > 
> > >  On Thu, 5 Nov 2015 at 01:10 Roger 

Re: [R-sig-Geo] How to perform cross-year date operations on rasters?

[Dominik Schneider]

Could you use water years for your dates? so Oct-Sep is the same year and
then subset based on month and water year.


On Thu, Nov 5, 2015 at 11:17 AM, Thiago V. dos Santos <
thi_vel...@yahoo.com.br> wrote:

> Thanks for your input Michael. With slight modifications on your
> suggestion, I almost got there:
>
> library(raster)
> library(zoo)
>
> # Create a rasterStack similar to cmip5 - same dimensions and layer names
> r <- raster(ncol=180, nrow=90)
> s <- stack(lapply(1:1825, function(x) setValues(r, runif(ncell(r)
> idx <- seq(as.Date("2010/1/1"), by = "day", length.out = 1825)
> s <- setZ(s, idx)
>
> # Separate layers with months of interest
> ldates <- format(getZ(s), "%m") %in% c("10", "11", "12", "01")
> s2 <- subset(s, which(ldates))
>
> # Apply function
> s3 <- zApply(s2, by=as.yearmon, fun = sum)
>
>
> The problem now is that the "isolated" january in the first year is also
> taken into account, and the "isolated" oct-nov-dec in the last year as well.
>
> Ideally the function would run only on the contiguous period:
> oct-nov-dec-jan, and not when at least one month is missing (for example in
> the first and last years of the series).
>
> Still possible to achieve this?
> Greetings,
>  -- Thiago V. dos Santos
>
> PhD student
> Land and Atmospheric Science
> University of Minnesota
>
>
>
> On Thursday, November 5, 2015 12:15 AM, Michael Sumner 
> wrote:
>
>
>
>
>
>
>
> On Thu, 5 Nov 2015 at 09:38 Thiago V. dos Santos 
> wrote:
>
> Dear all,
> >
> >Consider that I have a raster stack with daily values. How can I perform
> date operations covering a time interval that crosses years?
> >
> >For example, I want to sum the values from every october-to-january
> period in this sample raster:
> >
> >library(raster)
> >
> ># Create a rasterStack similar to cmip5 - same dimensions and layer names
> >r <- raster(ncol=180, nrow=90)
> >s <- stack(lapply(1:1825, function(x) setValues(r, runif(ncell(r)
> >
> ># Apply time stamps to raster
> >#x <- as.Date(c("2010-01-01","2014-12-31"),format="%Y-%m-%d")
> >#difftime(x[2], x[1], units="days")
> >idx <- seq(as.Date("2010/1/1"), by = "day", length.out = 1825)
> >s <- setZ(s, idx)
> >s
> >
> >
> >
>
>
> You can subset on the dates, by running a test on the date values:
>
> ldates <- format(getZ(s), "%m") %in% c("10", "11", "01")
>
> subsetting the object
>
> subset(s, which(ldates))
>
> and finally calculating what you want
>
> calc(subset(s, which(ldates)), sum)
>
> HTH
>
>
>
>
>  Thanks in advance,
> > -- Thiago V. dos Santos
> >
> >PhD student
> >Land and Atmospheric Science
> >University of Minnesota
> >
> >___
> >R-sig-Geo mailing list
> >R-sig-Geo@r-project.org
> >https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> ___
> 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] How to perform cross-year date operations on rasters?

[Thiago V. dos Santos]

Thanks for your input Michael. With slight modifications on your suggestion, I 
almost got there:

library(raster)
library(zoo)

# Create a rasterStack similar to cmip5 - same dimensions and layer names
r <- raster(ncol=180, nrow=90)
s <- stack(lapply(1:1825, function(x) setValues(r, runif(ncell(r)
idx <- seq(as.Date("2010/1/1"), by = "day", length.out = 1825)
s <- setZ(s, idx)

# Separate layers with months of interest
ldates <- format(getZ(s), "%m") %in% c("10", "11", "12", "01")
s2 <- subset(s, which(ldates))

# Apply function
s3 <- zApply(s2, by=as.yearmon, fun = sum)


The problem now is that the "isolated" january in the first year is also taken 
into account, and the "isolated" oct-nov-dec in the last year as well. 

Ideally the function would run only on the contiguous period: oct-nov-dec-jan, 
and not when at least one month is missing (for example in the first and last 
years of the series).

Still possible to achieve this? 
Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota



On Thursday, November 5, 2015 12:15 AM, Michael Sumner  
wrote:







On Thu, 5 Nov 2015 at 09:38 Thiago V. dos Santos  
wrote:

Dear all,
>
>Consider that I have a raster stack with daily values. How can I perform date 
>operations covering a time interval that crosses years?
>
>For example, I want to sum the values from every october-to-january period in 
>this sample raster:
>
>library(raster)
>
># Create a rasterStack similar to cmip5 - same dimensions and layer names
>r <- raster(ncol=180, nrow=90)
>s <- stack(lapply(1:1825, function(x) setValues(r, runif(ncell(r)
>
># Apply time stamps to raster
>#x <- as.Date(c("2010-01-01","2014-12-31"),format="%Y-%m-%d")
>#difftime(x[2], x[1], units="days")
>idx <- seq(as.Date("2010/1/1"), by = "day", length.out = 1825)
>s <- setZ(s, idx)
>s
>
>
>


You can subset on the dates, by running a test on the date values: 

ldates <- format(getZ(s), "%m") %in% c("10", "11", "01")

subsetting the object

subset(s, which(ldates))

and finally calculating what you want

calc(subset(s, which(ldates)), sum)

HTH



 
 Thanks in advance,
> -- Thiago V. dos Santos
>
>PhD student
>Land and Atmospheric Science
>University of Minnesota
>
>___
>R-sig-Geo mailing list
>R-sig-Geo@r-project.org
>https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] gUnaryUnion Not Dissolving Correctly

[Matt Strimas-Mackey]

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 

[R-sig-Geo] clusterR + rasterToPoints

[Guillermo E. Ponce-Campos]

Is it possible to take advantage of multi-cores with function
rasterToPoints for converting a stack of several rasters into a x,y,z data
frame ?

Thanks,
Guillermo

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo