An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20110310/d25ff335/attachment.pl>
CRS problem with Polyset2SpatialPolygons (maptools, PBSmapping)
2 messages · Delphine Antoniucci, Roger Bivand
On Thu, 10 Mar 2011, Delphine Antoniucci wrote:
Dear all, Im very new to using the maptools, PBSmapping and sp libraries in R, so Im sorry if this is something obvious for most of you. I have had trouble with passing from a PolySet object to SpatialPolygonsDataFrame object. Basically I got two SpatialPolygonsDataFrame (SP) objects. I need information about their intersections. So I used the function 'joinPolys' (PBSmapping) with the operation "INT" after having transformed my SP into Polysets. The problem is that I got an error message when I try to transform again the result into a SP, due to the CRS to be used.
No, it isn't obvious. PBSmapping "knows" unprojected (geographical coordinates) and UTM (projected coordinates) with a zone number. It does not provide facilities to store any other values.
Here is the syntax I have used:
aaa <- joinPolys(SpatialPolygons2PolySet(mapproj),
SpatialPolygons2PolySet(piece), operation="INT")
So if you examine the projection of:
attr(SpatialPolygons2PolySet(mapproj), "projection")
you'll see that it is set to "1", to signal that it is not "LL" or "UTM".
This is copied in joinPolys() to the output object.
PolySet2SpatialPolygons() was written assuming that the limitations of the
PolySet representation of coordinate reference systems should be
respected. For the future, I'll have it insert as.character(NA) and issue
a warning.
For now, I suggest using spTransform() in rgdal on the input objects to
transform them to CRS("+proj=longlat +ellps=GRS80") and back out on
return. I hope that no datum shift will take place, you'd need to check
whether the initial and final coordinates matched.
Hope this helps,
Roger
bbb <- as.PolySet(aaa) bbb <- PolySet2SpatialPolygons(bbb, close_polys=TRUE)
Erreur dans PolySet2SpatialPolygons(bbb, close_polys = TRUE) : unknown coordinate reference system I think this is due to a wrong projection attribute. By default, the r ratio is 1 (1 x unit per y unit)
str(aaa,max.level=2)
Classes ?PolySet? and 'data.frame': 14 obs. of 5 variables:
$ PID: int 4651 4651 4651 4651 4651 4651 4651 4651 4651 4651 ...
$ SID: int 1 1 1 1 1 1 1 1 1 1 ...
$ POS: int 1 2 3 4 5 6 7 8 9 10 ...
$ X : num 485038 485085 485080 485084 485058 ...
$ Y : num 6506613 6506550 6506531 6506529 6506441 ...
- attr(*, "projection")= chr "1"
But the original CRS of the SPs ('mapproj' and 'piece') is:
print(proj4string(mapproj))
[1] " +proj=lcc +lat_1=44 +lat_2=49 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs" Would you have any suggestion to stay in this CRS? Should I spTransform() my SP in another CSR that is known by PolySet2SpatialPolygons? Or would you recommend an easier way of getting the information Im looking for without using the class PolySet? Any help/hint will be very appreciated! Many thanks, Delphine [[alternative HTML version deleted]]
Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no