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?
library(spatstat)
library(maptools)
library(lattice)
library(sp)
library(foreign)
require("rgdal")
div0A <- readOGR(dsn="Archive", layer="Projections")
> summary(div0A)
Object of class SpatialPolygonsDataFrame
Coordinates:
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.
Answer
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.
No comments:
Post a Comment