I am trying to figure out how I can transform my unit of length from decimal degree to km when I import the following simple ArcGIS polygon shape. Ultimately, I would like to have this information available when I assign this spatial polygon to my window in spatstat [i.e., Window area = xx square unit & Unit of length: xx km].

I did came across this post: R-sig-geo

This is what I've done so far, but I don't understand why/if I need to do so. Doesn't my div0A SpatialPolygonsDataFrame already has my CRS argument defined? Why do I need to spTransform then?

div0A <- readOGR(dsn="Archive", layer="Projections")
> summary(div0A)
Object of class SpatialPolygonsDataFrame
     min       max
x -80.00 -57.63430
y  66.25  78.16666
Is projected: FALSE
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Data attributes:
       Id      NPAzimutha         UTM20            UTM19           AlberEA      
 Min.   :0   Min.   :348877   Min.   :349232   Min.   :349162   Min.   :348656  
 1st Qu.:0   1st Qu.:348877   1st Qu.:349232   1st Qu.:349162   1st Qu.:348656  
 Median :0   Median :348877   Median :349232   Median :349162   Median :348656  
 Mean   :0   Mean   :348877   Mean   :349232   Mean   :349162   Mean   :348656  
 3rd Qu.:0   3rd Qu.:348877   3rd Qu.:349232   3rd Qu.:349162   3rd Qu.:348656  
 Max.   :0   Max.   :348877   Max.   :349232   Max.   :349162   Max.   :348656 

proj4string(div0A) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
xxlcc <- spTransform(div0A, CRS("+init=epsg:4326"))
res1 <- lapply(slot(xxlcc, "polygons"), function(x)
  sapply(slot(x, "Polygons"), slot, "area"))

> summary(unlist(res1))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  89.72   89.72   89.72   89.72   89.72   89.72

newW<-as(xxlcc, "owin")
> summary(newW)
Window: polygonal boundary
single connected closed polygon with 29 vertices
enclosing rectangle: [-80, -57.6343]x[66.25, 78.16666]units
Window area =  89.7204 square units

Does this look correct? How can I specify/know my unit of length in km? Moreover, this window is ~Baffin Bay, west of Greenland and I am wondering if I should be using a different datum? E.g. EPSG Projection 4747 - GR96? Although, I will have to create ppp object, for which the coordinates have been taken with GPS satellite navigation system.

I would appreciate if someone can walk me through this or suggest an online example or further readings. I read chapter 4: Spatial Data Import and Export in Applied Spatial Data Analysis with R (2008), but I'm obviously very new to all of this.

  • 1
    Welcome to gis.stackexchange. Please don't put tags into the thread title. Thanks! – underdark Apr 4 '12 at 14:28
  • 1
    Have you looked at geosphere? – whuber Apr 4 '12 at 16:05

EPSG 4326 is a geographic coordinate system and as such, has decimal degrees as units of length. If you want meters, you'd need a projected coordinate system. That's why you need to do an spTransform. To convert your coordinate system from a geographic to a projected one.

You already had the right commands, if I'm not mistaken. You just need the right projection. Baffin Island should be around EPSG 2960. You might want to check it out. Should work fine once you plug it into spTransform.

Your Answer

