[R-sig-Geo] intersection of objects SpatialPolygonsDataFrame

Agustin Lobo Agustin.Lobo at ija.csic.es
Wed Apr 9 20:26:24 CEST 2008


Thanks,

I get a problem in the inverse transform:
 > delme <- 
PolySet2SpatialPolygons(SpatialPolygons2PolySet(absUTMpolysHABS2))
Error in 
PolySet2SpatialPolygons(SpatialPolygons2PolySet(absUTMpolysHABS2)) :
   unknown coordinate reference system

probably because
 > absUTMpolysHABS2 at proj4string
CRS arguments: NA

despite the fact that I do have a *.prj file in the directory
where I import with  rgdal (thought that rgdal would read the *.prj file
if present).

I'm working on this (proj4 strings are not straightforward for ED50). 
Meanwhile wanted to run the example in the PolySet2SpatialPolygons
help page but:

 > nor_coast_poly_sp <- map2SpatialPolygons(nor_coast_poly, IDs=IDs,
+  proj4string=CRS("+proj=longlat +datum=wgs84"))
Error in CRS("+proj=longlat +datum=wgs84") :
   unknown elliptical parameter name

does anybody know the right proj4string in that case?

Agus


Andrew Niccolai escribió:
> Try 
> 
> ?PolySet2SpatialPolygons
> 
> This function will bring the newly intersected polygon back to class SPDF
> and then you should be able to attach the corresponding data.frame to the
> new polygon object.
> 
> Andrew Niccolai
> Doctoral Candidate
> Yale School of Forestry
> 
> 
>  
> -----Original Message-----
> From: r-sig-geo-bounces at stat.math.ethz.ch
> [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf Of Agustin Lobo
> Sent: Wednesday, April 09, 2008 9:48 AM
> To: r-sig-geo at stat.math.ethz.ch
> Subject: [R-sig-Geo] intersection of objects SpatialPolygonsDataFrame
> 
> Dear list,
> 
> I have two objects of class "SpatialPolygonsDataFrame"
> and I want to get the intersection. The first object (x)
> are circles (approached as polygons) with a fake data
> frame. The circles could be just in a SpatialPolygons
> object, I just put the fake data frame to get the same
> class as I would have got from importing a shp file via readOGR.
> The second object (absUTMpolysHABS2) is a map of habitats
> (imported from a shp file
> via readOGR). I want the circles to intersect the map
> so that I can calculate the area of each habitat in
> each circle. I'm not sure if I'm on the right path:
> 
> It seemed to me that I had to use function joinPolys()
> with operation="INT" from package PBSmapping,
> 
> a <- (joinPolys(SpatialPolygons2PolySet(x), 
> SpatialPolygons2PolySet(absUTMpolysHABS2)))
> 
> a looks good (circles subdivided into polygons)
> but SpatialPolygons2PolySet() makes me loss the information
> in the data frame of the map, so I do not know the habitat
> for each polygon within each circle.
> 
> Any way to relate the new polygons to the data frame of the map?
> Or should I use an entirely different apporach in R?
> Or should I export x to a GIS and do this operation there?
> 
> Thanks!
> 
> Agus
> 

-- 
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
LLuis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
email: Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster




More information about the R-sig-Geo mailing list