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

frankma md.franklin at gmail.com
Mon Aug 25 11:50:40 CEST 2014


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/

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 
library(rgdal)      # for map projection work; also loads sp
library(PBSmapping) # for GIS_like geospatial object manipulation / anslysis
iIRLluding poly
# 
# 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"))

#
# 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"))
#
# 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"
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.*

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.



More information about the R-sig-Geo mailing list