[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