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.

Here is the syntax I have used:

>   aaa <- joinPolys(SpatialPolygons2PolySet(mapproj),
SpatialPolygons2PolySet(piece), operation="INT")

>   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]]

