Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

So, the usual way we read a shapefile in R is via the maptools package, like this:

sfdata <- readShapeSpatial("/path/to/my/shapefile.shp", proj4string=CRS("+proj=longlat"))

However, I have a use case whereby I don't have a shapefile.shp but instead I have a series of polygon coordinates

16.484375 59.736328125,17.4951171875 55.1220703125,24.74609375 55.0341796875,22.5927734375 61.142578125,16.484375 59.736328125

and its corresponding projection

coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

How do I "instantiate" sfdata (which will be a "polygon object") directly from this data? (without going in a roundabout way of creating a shapefile with these data and then reading from the newly created shapefile)

share|improve this question

1 Answer 1

up vote 14 down vote accepted

First get the coordinates into a 2-column matrix:

> xym
         [,1]     [,2]
[1,] 16.48438 59.73633
[2,] 17.49512 55.12207
[3,] 24.74609 55.03418
[4,] 22.59277 61.14258
[5,] 16.48438 59.73633

Then create a Polygon, wrap that into a Polygons object, then wrap that into a SpatialPolygons object:

> library(sp)
> p = Polygon(xym)
> ps = Polygons(list(p),1)
> sps = SpatialPolygons(list(ps))

The reason for this level of complexity is that a Polygon is a simple ring, a Polygons object can be several rings with an ID (here set to 1) (so is like a single feature in a GIS) and a SpatialPolygons can have a CRS. Ooh, I should probably set it:

> proj4string(sps) = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

If you want to turn it into a SpatialPolygonsDataFrame (which is what comes our of readShapeSpatial when the shapefile is polygons) then do:

> data = data.frame(f=99.9)
> spdf = SpatialPolygonsDataFrame(sps,data)
> spdf

giving this:

> summary(spdf)
Object of class SpatialPolygonsDataFrame
Coordinates:
       min      max
x 16.48438 24.74609
y 55.03418 61.14258
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Data attributes:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   99.9    99.9    99.9    99.9    99.9    99.9 
share|improve this answer
    
+1 Very nice, clear exposition. It's great to see the code broken up by explanations rather than being offered as a monolithic block! –  whuber Dec 30 '11 at 14:56
    
Excellent ... great to see how these objects are put together! Need to see more of the R help pages written clearly like this. –  Simbamangu Dec 30 '11 at 17:26
    
Its something I have to re-teach myself every time I want to do it, so I take any opportunity to teach other people! –  Spacedman Dec 31 '11 at 13:13
    
excellent... how would I go about adding multiple unique id (f) polygons to the data frame? –  mga Jun 26 '13 at 15:39
    
@Spacedman: can you please explain, how you convert the polygon coordinates in the question to the xym matrix in your answer? –  Stefan Aug 19 '13 at 10:51

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.