Coordinates to GIS with R
En muchos casos, una capa de puntos viene dada en forma de excel o archivo .csv donde 2 columnas definen las coordenadas de estos. Si estas fuentes son enormes (> 100 millones de líneas) la cosa puede ponerse chunca en QGIS (los esriadictos quizá no tengan ese problema, como siempre). R tiene la solución, como siempre.
En este caso no hay mucho más que dar respuesta a esta pregunta de GIS StackExchange: What’s the fastest way to convert large csvs to point shapefile
En este caso, se proponen dos métodos (1) transformar todos los archivo .csv localizados en un directorio y enviarlos con epsg a una ubicación y (2) transformar un enorme .csv (en este caso si imita uno multiplicando un data.frame) en una archivo gpkg. Estos archivos tienen ciertas ventajas respecto a los shp, incluída la mayor cantidad de líneas escribibles en un mismo archivo.
Se debe considerar que el proceso parte de una función básica de la librería landtools, de la cual solo será necesaria la función dftocoords. En esta se transforma una tabla de coordenadas en un objeto espacial, permitiendo además la transformación entre los epsg de input y output (por ejemplo, si se tuviera un csv con coordenadas geográficas en grados y se quisiera aplicar una transformación a UTM).
El proceso en R
library(landtools)
library(mapview)
library(stringr)
library(sf)
#------------------------------------------------------------------------
# METHOD 1 --> CSV'S TO SHP
#------------------------------------------------------------------------
# List cvs files
csvs <- list.files("./data/", '.csv', full.names = TRUE)
# Read all csv files
dfs <- lapply(csvs, function(x) read.csv(x, sep = ';'))
# See where your coords are stored (next AtriCoords parameter)
names(dfs[[1]])
# TRANSFORM CSV (data.frame) INTO SPATIAL
layers <- lapply(dfs, function(x) dftopoints(x,
"EPSG:25829",
"EPSG:25829",
AtriCoords = c("x", "y")))
# Export the layers, setting the name as the original name
outnames <- paste0(str_replace(basename(csvs), '.csv', ''), '.shp')
outfiles <- paste0('./output/', outnames)
for(i in 1:length(layers)){st_write(layers[[i]],
outfiles[[i]], append = FALSE)}
# Check the results
# mapview(layers[[1]], zcol = 'Altura..m.')
#------------------------------------------------------------------------
# METHOD 2 --> NOT PARALLEL PROCESS
#------------------------------------------------------------------------
#------------------------------------------------------------------------
# Create a huge (data.frame) to imitate a huge csv
x <- dfs[[1]]
while(dim(x)[[1]] < 100000000){x <- rbind(x, x)}
str(x)
# export x for testing
# write.csv(x, file = './data/HUGECSV.csv')
# x <- read.csv('./data/HUGECSV.csv')
# TRANSFORM CSV (data.frame) INTO SPATIAL
system.time({
lay <- dftopoints(x,
"EPSG:25829",
"EPSG:25829",
AtriCoords = c("x", "y"))
})
# Export file
system.time({
fname <- paste0("./output/", 'try01.pgkg')
st_write(lay, fname, append = FALSE)
})
