[R-sig-Geo] Error message when dissolving polygons

Mon Aug 25 14:39:06 CEST 2014

In your subsequent message, you claim that by setting LL to UTM, you 
resolve the problem. This is wrong, as I'll show below. You have not said 
what you mean by "dissolve" either, the help pages of both 
maptools::unionSpatialPolygons and equivalently rgeos::gUnaryUnion show 
exactly what to do, but I don't think that this is what you actually mean.

On Mon, 25 Aug 2014, frankma wrote:

> My apologies Roger, thank you for letting me know the mailing list etiquette. 
>
> I understand the gist of what you are saying but I'm still unsure which
> projections I am mixing up. I am using the shapefile available here:
> http://www.cso.ie/en/census/census2011boundaryfiles/

The files at this source are:

"The shapefiles are available in the Irish Grid Reference System (TM65 / 
Irish Grid - EPSG Projection 29902)"

so are projected, but not UTM.

>
> My code and error are below. I realise that it may seem obvious to you but
> it isn't to me, I've spent hours trying to "dissolve polygons" and there is
> something I am not understanding here, I am sorry for this. Appreciate your
> response so far and any further help is very much appreciated, thank you. In
> terms of affiliations I am currently trying to complete an MSc dissertation
> at the University of Sheffield (UK) - I also work full time.
>
>
> library(maptools)   # for geospatial services; also loads foreign and sp
> library(gpclib)     # General Polygon Clipping library

Neither of maptools nor gpclib are needed. The gpclib package should 
never, ever be used - use polyclip instead, or rgeos if your objects 
inherit sp classes.

> library(rgdal)      # for map projection work; also loads sp

Good!

> library(PBSmapping) # for GIS_like geospatial object manipulation / anslysis
> iIRLluding poly
> #

Not needed.

> # Part 1:
> # First example, derived from the maptools unionSpatialPolygon() method
> documentation in the R 
> # maptools user manual: We calculate and sum the area of each polygon in the
> Ireland map. 
> # Then, use unionSpatialPolygons() to dissolve the county polygon boundries,
> and assign each
> # polygon's area to one of four regions, based on longitude thresholds.
> #
> print("start of demo. Hit key...")
> browser()
> print("reading and transforming Ireland Shapefile...")
> #
> IRLBase <- readShapePoly("Shape/counties.shp",proj4string=CRS("+proj=aea
> +ellps=GRS80 +datum=WGS84"))

Here you are setting the projection (coordinate reference system, CRS) to 
AEA when in all probability the input was as the data source says. All the 
proj4string= argument does is to insert the given CRS into the object; it 
says - I know that this is correct; it does not project, because that 
needs a from CRS and a to CRS.

>

The example in maptools uses readShapePoly() to avoid depending on rgdal, 
but only rgdal should be used in real work (rgdal reads the prj file in 
shapefiles to set the CRS directly); the maptools functions will probably 
be disabled shortly tooblige users to make the right choices:

> IRLBase <- readOGR(".", "Census2011_Admin_Counties_generalised20m")
OGR data source with driver: ESRI Shapefile
Source: ".", layer: "Census2011_Admin_Counties_generalised20m"
with 34 features and 20 fields
Feature type: wkbPolygon with 2 dimensions
> proj4string(IRLBase)
[1] "+proj=tmerc +lat_0=53.5 +lon_0=-8 +k=1.000035 +x_0=200000 +y_0=250000 
+datum=ire65 +units=m +no_defs +ellps=mod_airy 
+towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15"
> library(rgeos)
rgeos version: 0.3-6, (SVN revision 450)
  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
  Polygon checking: TRUE
> gA <- gArea(IRLBase, byid=TRUE)

gives the entity areas in square metres; note that the actual areas will 
depend on the degree of line generalisation.

With this file, which may not be the one you are using, you may also need 
to use the encoding="latin1" argument in readOGR() to get from D\xfan 
Laoghaire-Rathdown to Dún Laoghaire-Rathdown.

> #
> # Transform the polygons (which were read in as unprojected geographic
> coordinates)
> # to an Albers Equal Area projection.
> #
> IRLProj = spTransform(IRLBase,CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))

You've set the CRS to AEA on reading, so this is a no-op.

> #
> # Convert to a PolygonSet for compatability with PBSmapping package
> routines.
> #
> IRLProjPS = SpatialPolygons2PolySet(IRLProj)
> #
> plotPolys(IRLProjPS, proj =
> TRUE,col="wheat1",xlab="longitude",ylab="latitude")
> #
> # polygon area calculations
> #
> print("Calculating Ireland polygon areas...")
> attr(IRLProjPS, "projection") <- "LL"

But you set the CRS to AEA, which is neither UTM not LL - this is simply 
wrong.

> IRLPolyAreas = calcArea(IRLProjPS,rollup=1)
>
> *Error in calcArea(IRLProjPS, rollup = 1) :
>  To calculate the areas of polygons defined by longitude-latitude
> coordinates, this routine first projects them in UTM.
>
> Attempted to automatically calculate the missing 'zone' attribute, but
> that failed because the mean longitude falls outside the range
> -180 < x <= 180.*

As I said before, the coordinates in the data object are already 
projected. By your setting the PolySet to UTM, you get some numbers out 
(effectively the same as by doing it properly, because PBSmapping then 
calculates planar areas in the metric of the coordinates, rather than 
projecting on the fly from LL to UTM and then grapping planar areas). But 
it is certainly not UTM.

Which leaves the muddle about what you mean by dissolving; I can guess 
that you need to define the longitude thresholds as polygons in the same 
projection as the polygons, intersect the two using
rgeos::gIntersection(..., byid=TRUE), take rgeos::gArea(..., byid=TRUE) of 
the output, and pick apart the row.names() of the output object to see how 
many square metres of each county belongs to each longitude threshold.

Hope this clarifies,

Roger

>
> Thanks,
> Matt
>
>
>
> Roger Bivand wrote
>> On Mon, 25 Aug 2014, frankma wrote:
>> 
>>> Hi,
>>>
>>> I'm having the same issue. Did you manage to solve this?
>> 
>> Please always quote the problem verbatim, providing a reproducible 
>> example. Even if anyone cares to open nabble to find the link, nobody can 
>> actually know whether your problem is the same. The reason we ask for 
>> reproducible examples is that the effort needed by the questioner to 
>> create a reproducible example means that he/she has to look more closely 
>> at - here "I’m following a tutorial found in the web" - what they are 
>> trying to do, and very often see their mistake themselves. Learning to 
>> debug one's own work is not easy, takes time, but wastes less of other 
>> peoples' time than asking lazy questions - see the posting guide:
>> 
>> http://www.r-project.org/posting-guide.html
>> 
>> It is also very helpful to see an affiliation, so that anyone caring to 
>> answer knows where to start. This is especially important when the 
>> sender's address is non-institutional.
>> 
>> In the original question:
>> 
>> https://stat.ethz.ch/pipermail/r-sig-geo/2014-April/020764.html
>> 
>> the obvious error was to think that the PolySet class in PBSmapping 
>> understands other projections than UTM:
>> 
>> # library(maptools)
>> # zm.* not available
>>> input = readShapePoly("zm.shp", proj4string=CRS("+proj=lcc +ellps=GRS80 
>> +datum=WGS84"))
>>> projectedpolygons = spTransform(input, CRS("+proj=lcc +ellps=GRS80 
>> +datum=WGS84"))
>>> zmprojPS = SpatialPolygons2PolySet(projectedpolygons)
>>> attr(zmprojPS, "projection") <- "LL"
>> 
>> I assume that no answer was posted because the error message was totally 
>> obvious, that is that lcc is Lambert conformal conic, but LL is the code 
>> for geographical, unprojected, coordinates, so PBSmapping could not 
>> project to UTM to calculate areas as:
>> 
>> "the mean longitude falls outside the range -180 < x <= 180."
>> 
>> If you are also mixing up projections, then this also explains what you 
>> are seeing.
>> 
>> Hope this clarifies,
>> 
>> Roger
>> 
>>>
>>> Thanks,
>>> Matt
>>>
>>>
>>>
>>> --
>>> View this message in context:
>>> http://r-sig-geo.2731867.n2.nabble.com/Error-message-when-dissolving-polygons-tp7586109p7586979.html
>>> Sent from the R-sig-geo mailing list archive at Nabble.com.
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> 
>
>> R-sig-Geo@
>
>>> 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.Bivand@
>
>> 
>> _______________________________________________
>> R-sig-Geo mailing list
>
>> R-sig-Geo@
>
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> Roger Bivand wrote
>> On Mon, 25 Aug 2014, frankma wrote:
>> 
>>> Hi,
>>>
>>> I'm having the same issue. Did you manage to solve this?
>> 
>> Please always quote the problem verbatim, providing a reproducible 
>> example. Even if anyone cares to open nabble to find the link, nobody can 
>> actually know whether your problem is the same. The reason we ask for 
>> reproducible examples is that the effort needed by the questioner to 
>> create a reproducible example means that he/she has to look more closely 
>> at - here "I’m following a tutorial found in the web" - what they are 
>> trying to do, and very often see their mistake themselves. Learning to 
>> debug one's own work is not easy, takes time, but wastes less of other 
>> peoples' time than asking lazy questions - see the posting guide:
>> 
>> http://www.r-project.org/posting-guide.html
>> 
>> It is also very helpful to see an affiliation, so that anyone caring to 
>> answer knows where to start. This is especially important when the 
>> sender's address is non-institutional.
>> 
>> In the original question:
>> 
>> https://stat.ethz.ch/pipermail/r-sig-geo/2014-April/020764.html
>> 
>> the obvious error was to think that the PolySet class in PBSmapping 
>> understands other projections than UTM:
>> 
>> # library(maptools)
>> # zm.* not available
>>> input = readShapePoly("zm.shp", proj4string=CRS("+proj=lcc +ellps=GRS80 
>> +datum=WGS84"))
>>> projectedpolygons = spTransform(input, CRS("+proj=lcc +ellps=GRS80 
>> +datum=WGS84"))
>>> zmprojPS = SpatialPolygons2PolySet(projectedpolygons)
>>> attr(zmprojPS, "projection") <- "LL"
>> 
>> I assume that no answer was posted because the error message was totally 
>> obvious, that is that lcc is Lambert conformal conic, but LL is the code 
>> for geographical, unprojected, coordinates, so PBSmapping could not 
>> project to UTM to calculate areas as:
>> 
>> "the mean longitude falls outside the range -180 < x <= 180."
>> 
>> If you are also mixing up projections, then this also explains what you 
>> are seeing.
>> 
>> Hope this clarifies,
>> 
>> Roger
>> 
>>>
>>> Thanks,
>>> Matt
>>>
>>>
>>>
>>> --
>>> View this message in context:
>>> http://r-sig-geo.2731867.n2.nabble.com/Error-message-when-dissolving-polygons-tp7586109p7586979.html
>>> Sent from the R-sig-geo mailing list archive at Nabble.com.
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> 
>
>> R-sig-Geo@
>
>>> 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.Bivand@
>
>> 
>> _______________________________________________
>> R-sig-Geo mailing list
>
>> R-sig-Geo@
>
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
>
>
> --
> View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Error-message-when-dissolving-polygons-tp7586109p7586981.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at 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.Bivand at nhh.no