[R-sig-Geo] CRS problem with Polyset2SpatialPolygons (maptools, PBSmapping)

[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