Question

I am trying to plot a route in R. The .kml file comes from Marble, a virtual globe application. It looks like this:

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://earth.google.com/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2">
<Document>
    <visibility>1</visibility>
    <Folder>
        <name>Route Request</name>
        <visibility>1</visibility>
        <Placemark>
            <name>Bari</name>
            <visibility>1</visibility>
            <ExtendedData>
                <Data name="routingVisited">
                    <value>false</value>
                </Data>
            </ExtendedData>
            <Point>
                <coordinates>16.8511800000,41.1177300000,8.0000000000</coordinates>
            </Point>
        </Placemark>
        <Placemark>
            <name>Trani</name>
            <visibility>1</visibility>
            <ExtendedData>
                <Data name="routingVisited">
                    <value>false</value>
                </Data>
            </ExtendedData>
            <Point>
                <coordinates>16.4153700000,41.2727300000,10.0000000000</coordinates>
            </Point>
        </Placemark>
    </Folder>
        <Document>
        <name>47.1 km (OSRM)</name>
        <visibility>1</visibility>
        <Placemark>
            <name>Route</name>
            <visibility>1</visibility>
            <LineString>
                <coordinates>16.8512040000,41.1174370000 16.8508380000,41.1173690000 16.8504820000,41.1172040000 16.8494450000,41.1164640000 16.8489240000,41.1160100000 16.8488120000,41.1158720000 16.8488340000,41.1158380000 16.8513550000,41.1143390000 16.8521900000,41.1143730000 16.8524160000,41.1143590000 16.8526040000,41.1143640000 16.8526740000,41.1143290000 16.8526970000,41.1142450000 16.8528600000,41.1126580000
…
…

I can read the file with read.OGR:

library(rgdal)
> r1 <- "bari-trani.kml"
> kml1 <- readOGR(r1, layer = "Route Request") # start and end
OGR data source with driver: KML 
Source: "bari-trani.kml", layer: "Route Request"
with 2 features and 2 fields
Feature type: wkbPoint with 3 dimensions
> kml2 <- readOGR(r1, layer = "47.1 km (OSRM)") # route information
OGR data source with driver: KML 
Source: "bari-trani.kml", layer: "47.1 km (OSRM)"
with 24 features and 2 fields
Feature type: wkbLineString with 2 dimensions

I can also plot the route:

> plot(kml1)
> lines(kml2, col = "red")

Route

I have also learned how to create maps. I use the Openstreetmap package:

library(OpenStreetMap)
> topleft <- c(42.08428793301783, 11.9741182832513)
> bottomright <- c(37.49482899208904, 19.1591768770013)
> map <- openmap(topleft, bottomright, type = "osm-bw")
> plot(map)

South of Italy

However, I am unable to get the route onto the map. Plotting the route after the map does not show anything. Is there a way to do this?

Était-ce utile?

La solution

Your map object has the following projection:

R> map$tiles[[1]]$projection
CRS arguments:
 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
+k=1.0 +units=m +nadgrids=@null +no_defs 

and bounding box:

R> map$bbox
$p1
[1] 1332953 5173614
$p2
[1] 2132790 4508306

Most of OSM, including the main tiling system, uses a spherical Mercator projection. Note that KML by specification uses only a single projection, EPSG:4326 (aka WGS84).

You can use spTransform (from packages sp/rgdal) to reproject your data (untested):

library("sp")
library("rgdal")

crs <- CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
+k=1.0 +units=m +nadgrids=@null +no_defs")
spTransform(kml1, crs)

See also here.

Licencié sous: CC-BY-SA avec attribution
Non affilié à StackOverflow
scroll top