Creating x and y distance coordinates for R from a

问题:

I would like to use a .kml track file to make a set of x, y coordinates for use in R.

What I have right now is a GoogleEarth track, which I believe is a LineString. I have heard that the rgdal package is usually what people use, but it doesn't work on Mac versions of R. If possible, I'd like to do this on a Mac, where I do the rest of my analyses. If necessary, I can do the conversion on R64 with Windows, and then bring the coordinates to my Mac, but that seems...clunky.

The beginning of the .kml code looks like this:

    <?xml version="1.0" encoding="UTF-8"?>  
<kml xmlns="http://www.opengis.net/kml/2.2"  
xmlns:gx="http://www.google.com/kml/ext/2.2"  xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">
    <Document>  
    <name>Perimeter_Track.kml</name>
<Placemark>
    <name>ACTIVE LOG</name>
    <LineString>
        <coordinates>
            -157.80736808,21.4323436,20.324951171875  

I want to convert it into x , y coordinates in kilometers from a point in my map. The finished product will be a line outline of a body of water, with species abundance data overlaid on it.

I have tried a couple of methods already:
1. Converting the .kml file into a .csv and importing it to r using read.csv;
2. Importing coordinates using getKMLcoordinates in the maptools package.

The problem with (1) is that when I try to convert the .kml coords into csv, I get an error in the converter program (kmlcsv) that says it can't read the file (I'm not sure why- the error logs aren't available).

When I try (2), I get coordinates that are arranged weirdly.
spa<-getKMLcoordinates("Perimeter_Track.kml", ignoreAltitude=TRUE)
summary(spa) returns:

 Length Class  Mode   

[1,] 128 -none- numeric

[2,] 242 -none- numeric

[3,] 34 -none- numeric

[4,] 126 -none- numeric

I believe this is because the .kml file is actually four separate tracks, separated by small gaps (i.e., where they turned the GPS off for a short time, then started again). Do I need to import these all separately in order to get the whole shape? If so, how do I do this?

I would like, eventually, to get this shape on a grid that is x by y km, where the coordinates are in km instead of GPS coords. If anyone has any insight into how to do this, I would love to hear from you!
Thanks very much in advance.

回答1:

Even though a precompiled package isn't available, you can still install rgdal from its source on a mac like follows:

1. Install the "GDAL complete" framework from http://www.kyngchaos.com/software/frameworks.

2. Add the locations of the programs you just installed to your unix path. In mac terminal, do:

   PATH=/Library/Frameworks/GDAL.framework/unix/bin:/Library/Frameworks/PROJ.framework/unix/bin$PATH
   

3. Download the source for the rgdal package from CRAN at http://cran.r-project.org/web/packages/rgdal/index.html.

4. Open R and build/install the rgdal package. Note that we have to specify the locations to some of the stuff we just installed.

   install.packages('~/Downloads/rgdal_0.7-1.tar.gz', repos=NULL, type='source', configure.args=c('--with-proj-include=/Library/Frameworks/PROJ.framework/unix/include', '--with-proj-lib=/Library/Frameworks/PROJ.framework/unix/lib'))
   

This installs fine on my Mac OS X 10.6. Good luck!

---

So the basic idea with your data might be:

library(rgdal)
library(maptools)

# Load KML coordinates
coords = getKMLcoordinates('data.kml')
coords = SpatialPoints(coords, CRS('+proj=longlat'))

# Load US Maps (get from www.gadm.org)
load('USA_adm1.RData') 
hawaii = gadm[gadm$NAME_1 == 'Hawaii', ]

# Transform coordinates
hawaii.proj = spTransform(hawaii, CRS=CRS('+init=epsg:2784 +units=km'))
coords.proj = spTransform(coords, CRS=CRS('+init=epsg:2784 +units=km'))

# Plot
dev.new(width=4, height=4)
plot(hawaii.proj, axes=T, xlim=c(450,550), ylim=c(0,60))
points(coords.proj, pch=16, col='red')

Great place to live!

回答2:

Once you've read something into an "sp" class object (probably a SpatialLinesDataFrame here) using readOGR from rgdal, you can transform it to a cartesian system from lat-long with the spTransform function.

Which system to transform it into depends on where on the earth it is. There's a bunch of standard transforms that depend on longitude called 'UTM' zones (Universal Transverse Mercator). Simply look up the zone for your longitude, find the EPSG code, and fire up spTransform.

For the UK, there's a standard grid system called the Ordnance Survey Grid, which has an EPSG code of 27700. So to transform something in lat-long (EPSG:4326) to OSGB metres, I do:

mapOS = spTransform(mapLL,CRS=CRS("+init=epsg:27700"))

There's lots of examples in the help for spTransform.

Note this is all great only if your data are on a small part of the earth...