Lectura de archivos de información geográfica y sistemas de referencias
Daniela Ballari
Universidad de Cuenca
daniela.ballari@ucuenca.edu.ec
Marzo, 2015
En esta lección aprenderá a: - Cargar archivos shp de tipo punto, linea y polígono
- Cargar archivos de texto, tipo csv
- Cargar archivos raster
- Explorar y transformar sistemas de referencia
- Asignar coordenadas
- Guardar archivos de texto y shp
Para comenzar instale y cargue la libreria gdal y sp, en caso que no los haya instalado antes. Luego defina el directorio de trabajo. Por favor actualice el directorio de acuerdo a su preferencia.
#install.packages("rgdal")
#install.packages("sp")
#install.packages("raster")
#install.packages("rgeos")
library(sp)## Warning: package 'sp' was built under R version 3.1.3library(rgdal)## Warning: package 'rgdal' was built under R version 3.1.3## rgdal: version: 0.9-1, (SVN revision 518)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.1, released 2014/09/24
## Path to GDAL shared files: C:/Users/dani/Documents/R/win-library/3.1/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: C:/Users/dani/Documents/R/win-library/3.1/rgdal/projsetwd("C:/curso_R/")Luego, descargue los datos y descomprimalos en el directorio de trabajo.
shell.exec("https://drive.google.com/uc?export=download&id=0B0ea6usixQHFTktDSC04MHNJY2M")Explore las opciones gdal.
ogrDrivers() #formatos soportados## name write
## 1 AeronavFAA FALSE
## 2 ARCGEN FALSE
## 3 AVCBin FALSE
## 4 AVCE00 FALSE
## 5 BNA TRUE
## 6 CSV TRUE
## 7 DGN TRUE
## 8 DXF TRUE
## 9 EDIGEO FALSE
## 10 ESRI Shapefile TRUE
## 11 Geoconcept TRUE
## 12 GeoJSON TRUE
## 13 Geomedia FALSE
## 14 GeoRSS TRUE
## 15 GML TRUE
## 16 GMT TRUE
## 17 GPSBabel TRUE
## 18 GPSTrackMaker TRUE
## 19 GPX TRUE
## 20 HTF FALSE
## 21 Idrisi FALSE
## 22 KML TRUE
## 23 MapInfo File TRUE
## 24 Memory TRUE
## 25 MSSQLSpatial TRUE
## 26 ODBC TRUE
## 27 ODS TRUE
## 28 OpenAir FALSE
## 29 OpenFileGDB FALSE
## 30 PCIDSK TRUE
## 31 PDF TRUE
## 32 PDS FALSE
## 33 PGDump TRUE
## 34 PGeo FALSE
## 35 REC FALSE
## 36 S57 TRUE
## 37 SDTS FALSE
## 38 SEGUKOOA FALSE
## 39 SEGY FALSE
## 40 SUA FALSE
## 41 SVG FALSE
## 42 SXF FALSE
## 43 TIGER TRUE
## 44 UK .NTF FALSE
## 45 VRT FALSE
## 46 Walk FALSE
## 47 WAsP TRUE
## 48 XLSX TRUE
## 49 XPlane FALSE#list.files() #Listar todos los archivos
ogrListLayers(".") #archivos georreferenciados shp en directorio de trabajo## [1] "estaciones_32717" "nxprovincias"
## [3] "prov_zona_estudio" "rio_l"
## [5] "rio_torrente" "vias"
## [7] "AL_CENTOS_SALUD_2012" "AL_CENTROS_EDUCATIVOS_2013"
## [9] "AREA_INUNDACION" "buffer_c_salud"
## [11] "observaciones_georef"
## attr(,"driver")
## [1] "ESRI Shapefile"
## attr(,"nlayers")
## [1] 11ogrInfo(".","prov_zona_estudio") #información del archivo prov_zona_estudio## Source: ".", layer: "prov_zona_estudio"
## Driver: ESRI Shapefile number of rows 7
## Feature type: wkbPolygon with 2 dimensions
## Extent: (557012.9 9445216) - (979220 9841834)
## CRS: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## LDID: 87
## Number of fields: 11
## name type length typeName
## 1 DPA_PROVIN 4 80 String
## 2 DPA_DESPRO 4 80 String
## 3 DPA_VALOR 0 10 Integer
## 4 DPA_ANIO 4 80 String
## 5 REI_CODIGO 4 80 String
## 6 REN_CODIGO 4 80 String
## 7 PEE_CODIGO 4 80 String
## 8 preccount 2 24 Real
## 9 precsum 2 24 Real
## 10 precmean 2 24 Real
## 11 id 0 10 IntegerCargue una capa shp de tipo polígono
Para ello se utiliza la función readOGR, de la siguiente manera: nombre_variable <- readOGR(directorio donde “.” es el working directory, nombre del shp (sin extensión .shp) entre comillas).
zona_estudio<-readOGR(".","prov_zona_estudio") #fuente: editado de INEC## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "prov_zona_estudio"
## with 7 features and 11 fields
## Feature type: wkbPolygon with 2 dimensionsExplore el sistemas de referencia
Explore la clase de la variable zona_estudio y el sistema de referencia asignado (CRS).
class(zona_estudio)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"proj4string(zona_estudio)## [1] "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"Si el archivo no tuviera CRS asignado debe utilizarse la función proj4string para asignar el CRS adecuado. En cambio, si se quisiera cambiar el sistema de referencia (transormar) debe utilizarse la función spTransform.
utm17 <- "+proj=utm +zone=17 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs" #utm zona 17 sur con datum WGS84
wgs84 <- "+proj=longlat +ellps=WGS84" #coordenadas geograficas WGS84
proj4string(zona_estudio) <- CRS(utm17) #asignación. Dará un error dado que los datos ya disponen de CRS asignado## Warning in `proj4string<-`(`*tmp*`, value = <S4 object of class structure("CRS", package = "sp")>): A new CRS was assigned to an object with an existing CRS:
## +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## without reprojecting.
## For reprojection, use function spTransform in package rgdalzona_estudio_wgs84 <- spTransform(zona_estudio,CRS(wgs84))#transformaciónAcceda a datos y coordenadas de la variable zona_estudio que es una clase SpatialPolygonsDataFrame.
zona_estudio@data #datos## DPA_PROVIN DPA_DESPRO DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO
## 0 01 AZUAY 0 2012 05 01
## 1 03 CAÃAR 0 2012 05 01
## 2 06 CHIMBORAZO 0 2012 02 01
## 3 07 EL ORO 0 2012 03 02
## 4 11 LOJA 0 2012 05 01
## 5 14 MORONA SANTIAGO 0 2012 05 03
## 6 19 ZAMORA CHINCHIPE 0 2012 05 03
## PEE_CODIGO preccount precsum precmean id
## 0 593 9 1307.4945 145.2772 1
## 1 593 4 660.2622 165.0656 3
## 2 593 9 1529.6052 169.9561 6
## 3 593 8 903.4966 112.9371 7
## 4 593 13 2250.3139 173.1011 11
## 5 593 30 5979.2642 199.3088 14
## 6 593 15 3008.0197 200.5346 19coordinates(zona_estudio) #coordenadas de centros de polígonos## [,1] [,2]
## 0 704262.6 9667236
## 1 720735.2 9718927
## 2 753811.6 9781653
## 3 631931.2 9610850
## 4 649169.9 9547764
## 5 832894.6 9716260
## 6 732696.0 9538694names(zona_estudio) #nombres de campos## [1] "DPA_PROVIN" "DPA_DESPRO" "DPA_VALOR" "DPA_ANIO" "REI_CODIGO"
## [6] "REN_CODIGO" "PEE_CODIGO" "preccount" "precsum" "precmean"
## [11] "id"Cargue una capa shp de puntos
estaciones<-readOGR(".","estaciones_32717") #estaciones metereológicas## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "estaciones_32717"
## with 100 features and 14 fields
## Feature type: wkbPoint with 2 dimensionsproj4string(estaciones)## [1] "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"class(estaciones)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"names(estaciones)## [1] "C1" "CODIGO" "NOMBRE.DE" "TIPO" "ZONA.HIDRO"
## [6] "Lat" "Long" "ALTURA" "PROVINCIA" "INSTITUCIO"
## [11] "FECHA.DE.I" "FECHA.DE.L" "x" "y"summary(estaciones)## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## coords.x1 586178.3 826452.7
## coords.x2 9558491.8 9738060.1
## Is projected: TRUE
## proj4string :
## [+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0]
## Number of points: 100
## Data attributes:
## C1 CODIGO NOMBRE.DE TIPO
## 1MP : 6 M012 : 1 ARENALES-COLA DE SAN PABLO: 1 AP: 2
## 2MP : 8 M032 : 1 ARENILLAS : 1 AR: 3
## 2P :23 M040 : 1 AYAPAMBA : 1 CO:18
## 2P - 3P: 2 M045 : 1 BALSAS : 1 CP:10
## 3MP :22 M050 : 1 BIBLIAN : 1 PG:13
## 3P :39 M060 : 1 BIBLIAN INECEL : 1 PV:54
## (Other):94 (Other) :94
## ZONA.HIDRO Lat Long ALTURA
## Min. :170.0 Min. :-3.993 Min. :-80.22 Min. : 4.0
## 1st Qu.:190.0 1st Qu.:-3.684 1st Qu.:-79.79 1st Qu.: 601.5
## Median :200.0 Median :-3.331 Median :-79.32 Median :1632.5
## Mean :229.8 Mean :-3.288 Mean :-79.34 Mean :1603.9
## 3rd Qu.:280.0 3rd Qu.:-2.866 3rd Qu.:-78.89 3rd Qu.:2585.0
## Max. :280.0 Max. :-2.368 Max. :-78.06 Max. :3450.0
##
## PROVINCIA INSTITUCIO FECHA.DE.I FECHA.DE.L
## Min. : 1.00 D A C : 1 1973 / 2 / 1 : 6 1958 / 4 / 1 : 1
## 1st Qu.: 1.00 F A E : 2 1981 / 2 / 1 : 6 1963 / 11 / 23 : 1
## Median : 7.00 INAMHI :44 1972 / 1 / 1 : 4 1964 / 11 / 30 : 1
## Mean : 5.25 INECEL :23 1963 / 1 / 1 : 2 1966 / 2 / 2 : 1
## 3rd Qu.: 7.00 INERHI : 1 1963 / 5 / 16: 2 1974 / 5 / 16: 1
## Max. :19.00 INOCAR : 1 (Other) :77 (Other) : 4
## PREDESUR:28 NA's : 3 NA's :91
## x y
## Min. :586178 Min. :9558492
## 1st Qu.:634111 1st Qu.:9592672
## Median :686568 Median :9631718
## Mean :684052 Mean :9636422
## 3rd Qu.:734490 3rd Qu.:9682809
## Max. :826453 Max. :9738060
## Seleccione determinadas columnas del SpatialPointsDataFrame.
estaciones2<- estaciones[,c(2,3,6,7,8,9)]#opción 1: número de columna
estaciones2<- estaciones[,c("CODIGO","NOMBRE.DE","Lat","Long", "x","y")]#opción 1: nombre de columna
names(estaciones2)<- c("estacion","nombre", "lat", "long", "x", "y") #asignar nuevos nombres a las columnasVisualice las Provincias con las estaciones superpuestas.
l1 = list("sp.points", estaciones2, pch = 19)
spplot(zona_estudio, "DPA_PROVIN", sp.layout = list(l1))
Cargue un archivo de texto
Lea datos de observaciones de lluvia de las estaciones en formato csv y revise si existen valores faltantes. La función table le ayudará a crear una tabla de frecuencias.
observaciones <- read.table("observaciones_2003.csv")
class(observaciones)## [1] "data.frame"str(observaciones)## 'data.frame': 8030 obs. of 7 variables:
## $ FECHA : Factor w/ 365 levels "01/01/2003","01/02/2003",..: 1 13 25 37 49 61 73 85 97 109 ...
## $ VALOR : num 0 0 0 0 0 0 0 0 0 0 ...
## $ STATUS : Factor w/ 3 levels "","N","T": 1 1 1 1 1 1 1 1 1 1 ...
## $ ESTACION: Factor w/ 22 levels "M012","M032",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ year : int 2003 2003 2003 2003 2003 2003 2003 2003 2003 2003 ...
## $ day : int 1 2 3 4 5 6 7 8 9 10 ...
## $ month : int 1 1 1 1 1 1 1 1 1 1 ...table(complete.cases(observaciones$VALOR))##
## FALSE TRUE
## 1786 6244observaciones2 <- observaciones[complete.cases(observaciones$VALOR),]#excluir valores vacios
table(complete.cases(observaciones2$VALOR))##
## TRUE
## 6244Guarde la variable observaciones2 como archivo txt.
write.table(observaciones2, "observaciones_no_missing.txt")Asigne coordenadas a las observaciones
Proceda a unir los data.frame estaciones2 y observaciones2 con la función mergepara que las observaciones hereden las coordenadas de sus respectivas estaciones. Las columnas a unir son observaciones2$ESTACION y estaciones2$estacion.
observ_ref <- merge(observaciones2, estaciones2,by.x = "ESTACION", by.y = "estacion")
class(observ_ref)## [1] "data.frame"names(observ_ref)## [1] "ESTACION" "FECHA" "VALOR" "STATUS" "year"
## [6] "day" "month" "nombre" "lat" "long"
## [11] "x" "y" "coords.x1" "coords.x2"La variable obtenida es de tipo data.frame, es decir no está georreferenciada. Sin embargo, ahora el data.frame contiene una columna llamada x e y que se utilizarán para georreferenciar con la ayuda de la función coordinates. Luego, se asignará un CRS ya que el SpatialPointsDataFrame no tiene asignado ninguno.
coordinates(observ_ref) <- ~x+y
class(observ_ref)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"proj4string(observ_ref) ## [1] NAproj4string(observ_ref) <- CRS(utm17)Visualice los puntos observados.
l1 = list("sp.points", observ_ref, pch = 19, col="green")
spplot(zona_estudio, "DPA_PROVIN", col="grey", sp.layout = list(l1))
Finalmente, guarde la variable georreferenciada como archivo shp.
writeOGR(observ_ref, dsn = ".", layer = "observaciones_georef", driver = "ESRI Shapefile",overwrite_layer = TRUE)Cargue una capa shp de líneas
rios<-readOGR(".","rio_torrente")#fuente IGM 50k## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "rio_torrente"
## with 2093 features and 5 fields
## Feature type: wkbLineString with 2 dimensionsproj4string(rios)## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"rios <- spTransform(rios,CRS(utm17))#transformación
class(rios)## [1] "SpatialLinesDataFrame"
## attr(,"package")
## [1] "sp"names(rios)## [1] "f_code" "nombre" "hyc" "cat_hidro" "soc"Dibuje capa de líneas superpuesta a las provincias.
l1 = list("sp.lines", rios, col="blue")
spplot(zona_estudio, "DPA_PROVIN", col="grey", sp.layout = list(l1))
Para mejorar esta visualización vamos a recortar la capa de rios para solo visualizar los rios en la zona de estudio. Para ello necesitaremos la libreria rgeos, por lo que procederemos a instalarla y cargarla. Rgeos contiene funciones de operaciones vectoriales, como el caso de intersección que es la que usaremos.
library(rgeos)## Warning: package 'rgeos' was built under R version 3.1.3## rgeos version: 0.3-8, (SVN revision 460)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Polygon checking: TRUEcroped <- gIntersection(zona_estudio,rios, byid = TRUE) # esta función llevará un tiempo en ejecutarse, ten paciencia.
summary(croped)## Object of class SpatialLines
## Coordinates:
## min max
## x 557265.7 942295.5
## y 9446773.2 9841350.0
## Is projected: TRUE
## proj4string :
## [+proj=utm +zone=17 +south +ellps=WGS84 +datum=WGS84 +units=m
## +no_defs +towgs84=0,0,0]class(croped)## [1] "SpatialLines"
## attr(,"package")
## [1] "sp"l1 = list("sp.lines", croped, col="blue")
spplot(zona_estudio, "DPA_PROVIN", col="grey", sp.layout = list(l1))
Cargue una capa raster
Primero leeremos la capa rater con la libreria geotiff, obteniendose como resultado una variable de clase ‘SpatialGridDataFrame’.Luego se utilizarán distintas opciones de visualización.
land_cover <- readGDAL("LC_5min_global_2012.tif") #Mapa global de uso del suelo, periodo 2001-2012 fuente http://glcf.umd.edu/data/lc/ ## LC_5min_global_2012.tif has GDAL driver GTiff
## and has 1776 rows and 4320 columnsclass(land_cover)## [1] "SpatialGridDataFrame"
## attr(,"package")
## [1] "sp"summary(land_cover)## Object of class SpatialGridDataFrame
## Coordinates:
## min max
## x -180 180
## y -64 84
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Grid attributes:
## cellcentre.offset cellsize cells.dim
## x -179.95833 0.08333333 4320
## y -63.95833 0.08333333 1776
## Data attributes:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 2.588 3.000 16.000str(land_cover)## Formal class 'SpatialGridDataFrame' [package "sp"] with 4 slots
## ..@ data :'data.frame': 7672320 obs. of 1 variable:
## .. ..$ band1: int [1:7672320] 0 0 0 0 0 0 0 0 0 0 ...
## ..@ grid :Formal class 'GridTopology' [package "sp"] with 3 slots
## .. .. ..@ cellcentre.offset: Named num [1:2] -180 -64
## .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
## .. .. ..@ cellsize : num [1:2] 0.0833 0.0833
## .. .. ..@ cells.dim : int [1:2] 4320 1776
## ..@ bbox : num [1:2, 1:2] -180 -64 180 84
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "x" "y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"image(land_cover)
spplot(land_cover) # lento pero funciona
Otra forma de cargar una capa raster es con la libreria raster. Una vez cargada, transformaremos su sistema de referencia, recortaremos la capa a nuestra zona de estudio con la función crop() y gardaremos la variable como archivo geotiff. Finalmente la variable se transforma a la clase ‘SpatialGridDataFrame’.
library(raster)## Warning: package 'raster' was built under R version 3.1.3land_cover2 <- raster( "LC_5min_global_2012.tif")
class(land_cover2)## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"land_cover2## class : RasterLayer
## dimensions : 1776, 4320, 7672320 (nrow, ncol, ncell)
## resolution : 0.08333333, 0.08333333 (x, y)
## extent : -180, 180, -64, 84 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : D:\DATOS\ejemplos_R\LC_5min_global_2012.tif
## names : LC_5min_global_2012
## values : 0, 255 (min, max)res(land_cover2) #resolución del raster## [1] 0.08333333 0.08333333extent(land_cover2) #extensión del raster## class : Extent
## xmin : -180
## xmax : 180
## ymin : -64
## ymax : 84proj4string(land_cover2) ## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"projection(land_cover2) ## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"#land_cover2_utm17 <- spTransform(land_cover2,CRS(utm17))#transformación utilizando libreria rgdal no funciona! se debe buscar otra alternativa.
land_cover2_utm17 <- projectRaster(land_cover2,crs=utm17)#transformación utilizando libreria raster si funciona## Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, :
## 1704 projected point(s) not finiteland_cover2_utm17## class : RasterLayer
## dimensions : 683, 1224, 835992 (nrow, ncol, ncell)
## resolution : 28500, 31500 (x, y)
## extent : -17495494, 17388506, -770656.6, 20743843 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=17 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0
## data source : in memory
## names : LC_5min_global_2012
## values : -3.648203, 20.38352 (min, max)land_cover2_croped <- crop(land_cover2_utm17, zona_estudio, snap="out")
plot(land_cover2_croped)
plot(zona_estudio, add = TRUE)
land_cover2sp<- as(land_cover2_croped, 'SpatialGridDataFrame')#transformar a libreria sp
image(land_cover2sp)
plot(zona_estudio, add = TRUE)
writeRaster(land_cover2_croped, filename="land_cover2_croped.tif", format="GTiff", overwrite=TRUE)## class : RasterLayer
## dimensions : 13, 16, 208 (nrow, ncol, ncell)
## resolution : 28500, 31500 (x, y)
## extent : 545005.8, 1001006, 9435343, 9844843 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : D:\DATOS\ejemplos_R\land_cover2_croped.tif
## names : land_cover2_croped
## values : 0, 13.99926 (min, max)Ha llegado al final de esta lección!