Error message when dissolving polygons
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@
-- 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@
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@
-- 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@
-- 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