¿Leyendo el archivo KML en R?

42

Estoy trabajando con grandes archivos .kml (hasta 10 Gb) y necesito una forma eficiente de leerlos en R. Hasta ahora, los he convertido en archivos de forma a través de QGIS y luego de vuelta a R con readShapePoly y readOGR (este último , por cierto, es ~ 1000 más rápido que el anterior). Idealmente, me gustaría cortar la etapa intermedia de QGIS, ya que es engorrosa y lenta.

¿Cómo leer archivos .kml directamente?

Yo veo esto también se puede hacer con readOGR . Desafortunadamente, no puedo ver cómo implementar el ejemplo trabajado (después de una larga preparación del archivo .kml:) xx <- readOGR(paste(td, "cities.kml", sep="/"), "cities"). Parece que "ciudades" aquí es el nombre de los objetos espaciales.

Roger Bivand admite que "cómo se descubre este nombre no es obvio, ya que el controlador KML en OGR lo necesita para acceder al archivo. Una posibilidad es:

system(paste("ogrinfo", paste(td, "cities.kml", sep="/")), intern=TRUE)

"

Pero esto tampoco funciona para mí. Aquí hay un archivo de prueba .kml para probarlo. Con él en mi directorio de trabajo, readOGR("x.kml", "id")genera este mensaje de error:

Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = use_iconv) : 
  Cannot open layer . 

Y system(paste("ogrinfo", "x.kml"), intern=TRUE)genera:

[1] "Had to open data source read-only."   "INFO: Open of `x.kml'"               
[3] "      using driver `KML' successful." "1: x (3D Polygon)"  

, que simplemente no entiendo.

¿Sería getKMLcoordinates{maptools} una alternativa válida?

También he intentado esto:

tkml <- getKMLcoordinates(kmlfile="x.kml", ignoreAltitude=T)
head(tkml[[1]])
tkml <- SpatialPolygons(tkml, 
                        proj4string=CRS("+init=epsg:3857"))

Las coordenadas se generan correctamente, pero mi intento de convertirlas nuevamente en un objeto poligonal falló con el siguiente mensaje:

Error in SpatialPolygons(tkml, proj4string = CRS("+init=epsg:3857")) : 
  cannot get a slot ("area") from an object of type "double"
RobinLovelace
fuente
1
Puede obtener las capas en el kml utilizando la función ogrListLayers de rgdal.
Mario Becerra

Respuestas:

37

Para leer un KML con el controlador OGR, debe asignarle el nombre del archivo y el nombre de la capa.

El comentario de Roger es que el nombre de la capa está oculto en el archivo KML y, a menos que sepa cómo se creó el KML, no puede inferir el nombre de la capa a partir del nombre del archivo KML.

Mirando su ejemplo KML, puedo ver:

<?xml version="1.0" encoding="utf-8" ?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<Document><Folder><name>x</name>
<Schema name="x" id="x">

Lo que me dice que el nombre de la capa es x, no id, y así:

> foo = readOGR("/tmp/x.kml", "x")
OGR data source with driver: KML 
Source: "/tmp/x.kml", layer: "x"
with 1 features and 2 fields
Feature type: wkbPolygon with 2 dimensions

funciona bien

Ahora, se puede tratar de obtener el nombre mediante el análisis del KML como XML usando un analizador R XML, o puede quizá tratar de leerlo en I como un archivo de texto hasta que encuentre la etiqueta con su nombre.

El otro enfoque es ejecutar el programa ogrinfo de línea de comandos que escupe los nombres de capa de un archivo KML:

$ ogrinfo /tmp/x.kml 
Had to open data source read-only.
INFO: Open of `/tmp/x.kml'
      using driver `KML' successful.
1: x (Polygon)

aquí mostrando que hay una capa de polígono llamada x.

Hombre espacial
fuente
Gracias por su respuesta Espaciado: resolvió el problema de inmediato. ¡Es una explicación clara como esta que me hace amar el intercambio de fichas! Una pregunta de 'punto extra': ¿podría usar el mismo comando para leer un subconjunto de datos (por ejemplo, el primer millón de polígonos)? De lo contrario, buscará dividir los enormes kmls con un programa externo.
RobinLovelace
2
KML siendo XML no está realmente diseñado para acceso aleatorio. La solución real es poner sus datos espaciales en una base de datos espacial y tener algunos índices espaciales de velocidad. Echa un vistazo a PostGIS.
Spacedman
Bien, buen plan: le he dicho al cliente que PostGIS es el camino a seguir para esos grandes datos, y estoy convencido de que es la opción correcta para el tipo de cosas que quiere hacer. ¡Buena excusa para que lo aprenda correctamente!
RobinLovelace
También existe la extensión espacial de sqlite , una base de datos basada en archivos, que no requeriría que instale un servicio y requiere menos configuración que PostGIS.
Frank
curiosamente systemen I necesaria path.expanden ~para ogrinfoal trabajo, a pesar de que funcionó bien en el camino sin expandir en la línea de comandos (macOS; Sys.which('ogrinfo')y which ogrinfodevuelto los mismos caminos)
MichaelChirico
5

Si desea hacer la forma alternativa usando maptool, esto debería funcionar:

tkml <- getKMLcoordinates(kmlfile="yourkml.kml", ignoreAltitude=T)
#make polygon
p1 = Polygon(tkml)
#make Polygon class
p2 = Polygons(list(p1), ID = "drivetime")
#make spatial polygons class
p3= SpatialPolygons(list(p2),proj4string=CRS("+init=epsg:4326"))

La clave aquí es que debe seguir un par de pasos para crear una clase de polígono espacial.

Visto
fuente
hola @Seen, he probado tu enfoque pero parece que no funciona. Tengo un error: Error en Polígono (tkml): las coordenadas deben ser una matriz de dos columnas> cabeza (tkml) [[1]] [1] -87.88141 30.49800 y lo tengo como una lista ... ¿crees que está bien convertir? lista de coordenadas a la matriz? tahnks!
maycca
1

No sé si esto sigue siendo un problema para alguien más, pero estuve corriendo en círculos por un tiempo con esto. Lo que finalmente funcionó para mí está abajo. Utiliza el XMLpaquete para llegar al xmlValuenodo correcto. Tuve que establecer el layerparámetro de readOGRal nombre de una de las carpetas dentro del archivo kml. Cuando configuro el layerparámetro en el archivo kml, obtengo el mismo error que RobinLovelace describió anteriormente.

A continuación se muestran muchas líneas de código que solo muestran cómo ver los distintos niveles de nodo del documento kml. Creo que esto será ligeramente diferente dependiendo de la fuente del kml. Pero debería poder usar la misma lógica para determinar el valor de parámetro correcto.

Además, he creado una lista de archivos KML por lo que se podría hacer fácilmente en una función que se podría poner en un lapply- do.callpar. Esto podría extraer datos de una larga lista de archivos kml. O bien, muchas subcarpetas dentro de un solo archivo kml, ya que parece readOGRque no pueden manejar múltiples subcarpetas en un archivo kml.

library(rgdal); library(XML)

# SET WORKING DIRECTORY FIRST!!
dir <- getwd()

kmlfilelist <- list.files(dir, pattern =".kml$", full.names=TRUE, recursive=FALSE)

doc0 <- xmlTreeParse(kmlfilelist[2], useInternal = TRUE)
rootNode0 <- xmlRoot(doc0)
rootName0 <- xmlName(rootNode0)
element1Name0 <- names(rootNode0)

nodeNames <- names(rootNode0[1][[1]])

# entire rootNode - kml Document level
rootNode0[[1]]

# 1st element of rootNode - kml file name
rootNode0[[1]][[1]] 

# 2nd element of rootNode - kml Style Map 
rootNode0[[1]][[2]] 

# 3rd element of rootNode - Style
rootNode0[[1]][[3]]

# 4th element of rootNode - Style
rootNode0[[1]][[4]] 

# 5th element of rootNode - kml Folder with data in it.
rootNode0[[1]][[5]] 

# 5th element 1st subelement of rootNode - kml Folder name with data in it. 
#  What to set readOGR() layer parameter to.
rootNode0[[1]][[5]][[1]] 

kmlfoldername <- xmlValue(rootNode0[[1]][[5]][[1]]) # Folder name to set = layer.

readOGR(dsn=kmlfilelist[2], layer =  kmlfoldername)
Arcilla
fuente
0

No sé si debería haber modificado mi respuesta anterior. Quizás, pero eso cubre algunas cosas que no están en esta respuesta, así que decidí dejarlo.

De todos modos, el siguiente código funciona bien para mí. Busca todos los xmlNodes en el archivo kml que se llaman "Carpeta" y luego establece el layerparámetro de readOGReso xmlValue. Probado en el directorio de trabajo con alrededor de 6 archivos kml separados. La salida es una lista de objetos SpatialDataFrames importados. Cada SpatialDataFrame puede ser fácilmente subconjunto de la lista.

Todavía no aborda los archivos kml con múltiples nodos de carpeta. Pero esa característica podría agregarse fácilmente con otra applyfunción anidada .

library(rgdal); library(XML)

# SET WORKING DIRECTORY FIRST!!
dir <- getwd()

kmlfilelist <- list.files(dir, pattern =".kml$", full.names=TRUE, recursive=FALSE)

ImportKml <- function (kmlfile) {
  doc0 <- xmlTreeParse(kmlfile, useInternal = TRUE)
  rootNode0 <- xmlRoot(doc0)
  rootName0 <- xmlName(rootNode0)
  element1Name0 <- names(rootNode0)

  kmlNodeNames <- unname(names(rootNode0[1][[1]]))
  kmlFolderNodeNum <- which(kmlNodeNames == "Folder")
  kmlFolderNodeName <- xmlValue(rootNode0[[1]][[kmlFolderNodeNum]][[1]])

  kmlIn <- readOGR(dsn=kmlfile, layer = kmlFolderNodeName)
}
ImportedKmls <- lapply(kmlfilelist, ImportKml)
Arcilla
fuente