Lectura de archivos de geoinformación y sistemas de referencia espacial
Daniela Ballari
Julio, 2018
Maestría en Hidrología
Universidad de Cuenca
http://www.ucuenca.edu.ec/maestria-ecohidrologia/
Daniela Ballari (PhD)
daniela.ballari@ucuenca.edu.ec
http://www.researchgate.net/profile/Daniela_Ballari/publications
Johanna Orellana Alvear (MSc.)
johanna.orellana@ucuenca.edu.ec
Curso completo en: http://rpubs.com/Johanna_Orellana_Alvear/MHidro_indice_2018
En esta lección aprenderás 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, shp y raster
Temario
A- Librerías
B- Carga una capa shp de tipo polígono
C- Explora el sistemas de referencia espacial
D- Carga una capa shp de puntos
E- Carga una capa shp de l?neas
F- Carga un archivo de texto csv
H- Georreferenciar observaciones meteorológicas
I- Guarda los archivos creados
A- Librerías
Para comenzar carga las librerías gdal,sp, rastery rgeos y define el directorio de trabajo.
library(sp)
library(rgdal)
library(raster)
library(rgeos)Explora algunas de las opciones de la libreria rgdal.
setwd("~/Documents/R_Workspace/Maestria_Hidrologia_2018/Dani_clase_adicional/datos/SRS")
ogrDrivers() # Listado de formatos soportados## name long_name
## 1 ARCGEN Arc/Info Generate
## 2 AVCBin Arc/Info Binary Coverage
## 3 AVCE00 Arc/Info E00 (ASCII) Coverage
## 4 AeronavFAA Aeronav FAA
## 5 AmigoCloud AmigoCloud
## 6 BNA Atlas BNA
## 7 CSV Comma Separated Value (.csv)
## 8 CSW OGC CSW (Catalog Service for the Web)
## 9 Carto Carto
## 10 Cloudant Cloudant / CouchDB
## 11 CouchDB CouchDB / GeoCouch
## 12 DGN Microstation DGN
## 13 DXF AutoCAD DXF
## 14 EDIGEO French EDIGEO exchange format
## 15 ESRI Shapefile ESRI Shapefile
## 16 ElasticSearch Elastic Search
## 17 GFT Google Fusion Tables
## 18 GML Geography Markup Language (GML)
## 19 GPKG GeoPackage
## 20 GPSBabel GPSBabel
## 21 GPSTrackMaker GPSTrackMaker
## 22 GPX GPX
## 23 GeoJSON GeoJSON
## 24 GeoRSS GeoRSS
## 25 Geoconcept Geoconcept
## 26 HTF Hydrographic Transfer Vector
## 27 HTTP HTTP Fetching Wrapper
## 28 Idrisi Idrisi Vector (.vct)
## 29 JML OpenJUMP JML
## 30 KML Keyhole Markup Language (KML)
## 31 MapInfo File MapInfo File
## 32 Memory Memory
## 33 ODS Open Document/ LibreOffice / OpenOffice Spreadsheet
## 34 OGR_GMT GMT ASCII Vectors (.gmt)
## 35 OGR_PDS Planetary Data Systems TABLE
## 36 OGR_SDTS SDTS
## 37 OGR_VRT VRT - Virtual Datasource
## 38 OSM OpenStreetMap XML and PBF
## 39 OpenAir OpenAir
## 40 OpenFileGDB ESRI FileGDB
## 41 PCIDSK PCIDSK Database File
## 42 PDF Geospatial PDF
## 43 PGDUMP PostgreSQL SQL dump
## 44 PLSCENES Planet Labs Scenes API
## 45 REC EPIInfo .REC
## 46 S57 IHO S-57 (ENC)
## 47 SEGUKOOA SEG-P1 / UKOOA P1/90
## 48 SEGY SEG-Y
## 49 SQLite SQLite / Spatialite
## 50 SUA Tim Newport-Peace's Special Use Airspace Format
## 51 SVG Scalable Vector Graphics
## 52 SXF Storage and eXchange Format
## 53 Selafin Selafin
## 54 TIGER U.S. Census TIGER/Line
## 55 UK .NTF UK .NTF
## 56 VDV VDV-451/VDV-452/INTREST Data Format
## 57 VFK Czech Cadastral Exchange Data Format
## 58 WAsP WAsP .map format
## 59 WFS OGC WFS (Web Feature Service)
## 60 XLSX MS Office Open XML spreadsheet
## 61 XPlane X-Plane/Flightgear aeronautical data
## 62 netCDF Network Common Data Format
## write copy isVector
## 1 FALSE FALSE TRUE
## 2 FALSE FALSE TRUE
## 3 FALSE FALSE TRUE
## 4 FALSE FALSE TRUE
## 5 TRUE FALSE TRUE
## 6 TRUE FALSE TRUE
## 7 TRUE FALSE TRUE
## 8 FALSE FALSE TRUE
## 9 TRUE FALSE TRUE
## 10 TRUE FALSE TRUE
## 11 TRUE FALSE TRUE
## 12 TRUE FALSE TRUE
## 13 TRUE FALSE TRUE
## 14 FALSE FALSE TRUE
## 15 TRUE FALSE TRUE
## 16 TRUE FALSE TRUE
## 17 TRUE FALSE TRUE
## 18 TRUE FALSE TRUE
## 19 TRUE TRUE TRUE
## 20 TRUE FALSE TRUE
## 21 TRUE FALSE TRUE
## 22 TRUE FALSE TRUE
## 23 TRUE FALSE TRUE
## 24 TRUE FALSE TRUE
## 25 TRUE FALSE TRUE
## 26 FALSE FALSE TRUE
## 27 FALSE FALSE TRUE
## 28 FALSE FALSE TRUE
## 29 TRUE FALSE TRUE
## 30 TRUE FALSE TRUE
## 31 TRUE FALSE TRUE
## 32 TRUE FALSE TRUE
## 33 TRUE FALSE TRUE
## 34 TRUE FALSE TRUE
## 35 FALSE FALSE TRUE
## 36 FALSE FALSE TRUE
## 37 FALSE FALSE TRUE
## 38 FALSE FALSE TRUE
## 39 FALSE FALSE TRUE
## 40 FALSE FALSE TRUE
## 41 TRUE FALSE TRUE
## 42 TRUE TRUE TRUE
## 43 TRUE FALSE TRUE
## 44 FALSE FALSE TRUE
## 45 FALSE FALSE TRUE
## 46 TRUE FALSE TRUE
## 47 FALSE FALSE TRUE
## 48 FALSE FALSE TRUE
## 49 TRUE FALSE TRUE
## 50 FALSE FALSE TRUE
## 51 FALSE FALSE TRUE
## 52 FALSE FALSE TRUE
## 53 TRUE FALSE TRUE
## 54 TRUE FALSE TRUE
## 55 FALSE FALSE TRUE
## 56 TRUE FALSE TRUE
## 57 FALSE FALSE TRUE
## 58 TRUE FALSE TRUE
## 59 FALSE FALSE TRUE
## 60 TRUE FALSE TRUE
## 61 FALSE FALSE TRUE
## 62 TRUE TRUE TRUElist.files() # Listado de archivos en el directorio de trabajo## [1] "LC_5min_global_2012.tif" "elevacion.tif"
## [3] "elevacion.tif.aux.xml" "estaciones_32717.dbf"
## [5] "estaciones_32717.prj" "estaciones_32717.qpj"
## [7] "estaciones_32717.shp" "estaciones_32717.shx"
## [9] "observaciones_2003.csv" "prov_zona_estudio.dbf"
## [11] "prov_zona_estudio.prj" "prov_zona_estudio.shp"
## [13] "prov_zona_estudio.shx" "rio_32717_zona_estudio.dbf"
## [15] "rio_32717_zona_estudio.prj" "rio_32717_zona_estudio.qpj"
## [17] "rio_32717_zona_estudio.shp" "rio_32717_zona_estudio.shx"
## [19] "rio_torrente.dbf" "rio_torrente.prj"
## [21] "rio_torrente.shp" "rio_torrente.shx"ogrListLayers(".") #Archivos georreferenciados shp en el directorio de trabajo## [1] "estaciones_32717" "prov_zona_estudio"
## [3] "rio_32717_zona_estudio" "rio_torrente"
## attr(,"driver")
## [1] "ESRI Shapefile"
## attr(,"nlayers")
## [1] 4ogrInfo(".","prov_zona_estudio") #iInformación del archivo prov_zona_estudio.shp## 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 12 10 Integer64
## 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 12 10 Integer64B- Cargar una capa shp de tipo polígono
Para ello se utiliza la función readOGR de la siguiente manera: variableResultado <- readOGR(<directorio>, <nombre del shp>).
setwd("~/Documents/R_Workspace/Maestria_Hidrologia_2018/Dani_clase_adicional/datos/SRS")
zona_estudio<-readOGR(".","prov_zona_estudio") #fuente de datos: editado INEC## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "prov_zona_estudio"
## with 7 features
## It has 11 fields
## Integer64 fields read as strings: DPA_VALOR idC- Explora el sistemas de referencia espacial
Explora 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. Por otro lado, si se quisiera cambiar el sistema de referencia (transformar) debe utilizarse la función spTransform. A continuación se incluyen ejemplos de las dos operaciones.
#Primero se definen los sistemas utm17 y wgs84. Para ello hay dos opciones:
# Opción 1 #utm zona 17 sur con datum WGS84
utm17 <-"+init=epsg:32717"
# Opción 2 #utm zona 17 sur con datum WGS84
utm17 <- "+proj=utm +zone=17 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
# Opción 1 #coordenadas geograficas WGS84
wgs84 <-"+init=epsg:4326"
# Opción 2 #coordenadas geograficas WGS84
wgs84 <- "+proj=longlat +ellps=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 spTransformzona_estudio_wgs84 <- spTransform(zona_estudio,CRS(wgs84))#transformación
zona_estudio_wgs84## class : SpatialPolygonsDataFrame
## features : 7
## extent : -80.48635, -76.693, -5.016157, -1.429887 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84
## variables : 11
## names : DPA_PROVIN, DPA_DESPRO, DPA_VALOR, DPA_ANIO, REI_CODIGO, REN_CODIGO, PEE_CODIGO, preccount, precsum, precmean, id
## min values : 01, AZUAY, 0, 2012, 02, 01, 593, 4, 660.2622, 112.9371, 1
## max values : 19, ZAMORA CHINCHIPE, 0, 2012, 05, 03, 593, 30, 5979.2642, 200.5346, 7Visualiza los datos y coordenadas de la variable zona_estudio que es una clase SpatialPolygonsDataFrame.
zona_estudio@data #datos de tabla de atributos## DPA_PROVIN DPA_DESPRO DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO
## 0 01 AZUAY 0 2012 05 01
## 1 03 CA\303\221AR 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"D- Carga una capa shp de puntos
setwd("~/Documents/R_Workspace/Maestria_Hidrologia_2018/Dani_clase_adicional/datos/SRS")
estaciones<-readOGR(".","estaciones_32717") #estaciones meteorológicas## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "estaciones_32717"
## with 40 features
## It has 14 fields
## Integer64 fields read as strings: ZONA_HIDRO ALTURA PROVINCIAproj4string(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 731094.3
## coords.x2 9561300.3 9667787.9
## Is projected: TRUE
## proj4string :
## [+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0]
## Number of points: 40
## Data attributes:
## C1 CODIGO NOMBRE_DE TIPO ZONA_HIDRO
## 2P : 2 M012 : 1 ARENILLAS : 1 AP: 2 170: 1
## 2P - 3P: 2 M032 : 1 BALSAS : 1 CO: 9 180:14
## 3MP :15 M040 : 1 CHACRAS : 1 CP: 2 190: 5
## 3P :21 M142 : 1 COCHAPAMBA-QUINGEO: 1 PG: 2 200:12
## M179 : 1 CUMBE : 1 PV:25 210: 1
## M180 : 1 EL SALADO-PREDESUR: 1 280: 7
## (Other):34 (Other) :34
## Lat Long ALTURA PROVINCIA INSTITUCIO
## Min. :-3.968 Min. :-80.22 2525 : 2 1 :12 INAMHI :23
## 1st Qu.:-3.724 1st Qu.:-79.77 60 : 2 11: 4 INECEL : 4
## Median :-3.478 Median :-79.62 10 : 1 19: 1 INERHI : 1
## Mean :-3.492 Mean :-79.55 1040 : 1 7 :23 PREDESUR:12
## 3rd Qu.:-3.284 3rd Qu.:-79.25 1100 : 1
## Max. :-3.004 Max. :-78.92 1126 : 1
## (Other):32
## FECHA_DE_I FECHA_DE_L x y
## 1972 / 1 / 1 : 3 1958 / 4 / 1 : 1 Min. :586178 Min. :9561300
## 1981 / 2 / 1 : 3 1966 / 2 / 2 : 1 1st Qu.:636528 1st Qu.:9588266
## 1963 / 5 / 16: 2 1974 / 1 / 15: 1 Median :653549 Median :9615409
## 1973 / 2 / 1 : 2 1984 / 1 / 1 : 1 Mean :661292 Mean :9613894
## 1978 / 8 / 7 : 2 NA's :36 3rd Qu.:694320 3rd Qu.:9636876
## (Other) :25 Max. :731094 Max. :9667788
## NA's : 3Selecciona determinadas columnas del SpatialPointsDataFrame.
estaciones2<- estaciones[,c(2,3,6,7,8,9)]#número de columna
head(estaciones2)#Selecciona las 6 primeras filas## CODIGO NOMBRE_DE Lat Long ALTURA PROVINCIA
## 1 M012 LA CUCA -3.492222 -80.05694 53 7
## 2 M032 SANTA ISABEL INAMHI -3.274444 -79.31278 1550 1
## 3 M040 PASAJE -3.329722 -79.78194 40 7
## 4 M142 SARAGURO -3.620556 -79.23222 2525 11
## 5 M179 ARENILLAS -3.560278 -80.05611 60 7
## 6 M180 ZARUMA -3.696944 -79.61611 1100 7Visualiza las Provincias con las estaciones superpuestas.
l1 = list("sp.points", estaciones2, pch = 19)
spplot(zona_estudio, "DPA_PROVIN", sp.layout = list(l1))
E- Carga una capa shp de líneas
setwd("~/Documents/R_Workspace/Maestria_Hidrologia_2018/Dani_clase_adicional/datos/SRS")
rios<-readOGR(".","rio_torrente")#fuente IGM 50k## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "rio_torrente"
## with 2093 features
## It has 5 fieldsproj4string(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"Dibuja 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 ríos para solo visualizar los ríos en la zona de estudio. Para ello necesitaremos la libreria rgeos (funciones gLimeMerge y gIntersection). rgeos contiene funciones de operaciones vectoriales, como el caso de intersecci?n que es la que usaremos.
rios2 <- gLineMerge(rios)#Primero las lineas se unen con merge para que la siguiente funcion gIntersection sea mas rápida.
croped <- gIntersection(zona_estudio,rios2, byid = TRUE)
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))
F- Carga un archivo de texto csv
Lee los datos de observaciones de lluvia de las estaciones en formato csv y revisa si existen valores faltantes (missing values). La función table ayuda a crear una tabla de frecuencias.
setwd("~/Documents/R_Workspace/Maestria_Hidrologia_2018/Dani_clase_adicional/datos/SRS")
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 ...summary(observaciones)## FECHA VALOR STATUS ESTACION
## 01/01/2003: 22 Min. : 0.000 :7471 M012 : 365
## 01/02/2003: 22 1st Qu.: 0.000 N : 142 M032 : 365
## 01/03/2003: 22 Median : 0.000 T : 52 M040 : 365
## 01/04/2003: 22 Mean : 1.988 NA's: 365 M142 : 365
## 01/05/2003: 22 3rd Qu.: 1.000 M179 : 365
## 01/06/2003: 22 Max. :91.600 M180 : 365
## (Other) :7898 NA's :1786 (Other):5840
## year day month
## Min. :2003 Min. : 1.00 Min. : 1.000
## 1st Qu.:2003 1st Qu.: 8.00 1st Qu.: 4.000
## Median :2003 Median :16.00 Median : 7.000
## Mean :2003 Mean :15.72 Mean : 6.526
## 3rd Qu.:2003 3rd Qu.:23.00 3rd Qu.:10.000
## Max. :2003 Max. :31.00 Max. :12.000
## table(complete.cases(observaciones$VALOR))##
## FALSE TRUE
## 1786 6244observaciones2 <- observaciones[complete.cases(observaciones$VALOR),]#excluir valores vacios
table(complete.cases(observaciones2$VALOR))##
## TRUE
## 6244G- Carga un archivo raster
Esto lo podemos hacer con dos opciones.La primera con la función readGDAL. Como resultado se obtiene una variable de clase SpatialGridDataFrame. La segunda forma es con la libreria raster y se obtiene una variable de clase RasterLayer. Una vez cargada, transformaremos su sistema de referencia, y recortaremos la capa a nuestra zona de estudio con la función crop. Finalmente la variable se transforma a la clase SpatialGridDataFrame. Se utilizarán distintas opciones de visualización.
setwd("~/Documents/R_Workspace/Maestria_Hidrologia_2018/Dani_clase_adicional/datos/SRS")
####Opci?n 1 readGDAL
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:
## band1
## Min. : 0.000
## 1st Qu.: 0.000
## Median : 0.000
## Mean : 2.588
## 3rd Qu.: 3.000
## Max. :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)
####Opción 2
land_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 : /Users/johannaorellana/Documents/R_Workspace/Maestria_Hidrologia_2018/Dani_clase_adicional/datos/SRS/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 : 84plot(land_cover2)
proj4string(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 <- projectRaster(land_cover2,crs=utm17)#transformación SRS## 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")
image(land_cover2_croped)
plot(zona_estudio, colorkey=FALSE, add = TRUE)
land_cover2sp<- as(land_cover2_croped, 'SpatialGridDataFrame')#transformar a libreria sp
image(land_cover2sp)
plot(zona_estudio, add = TRUE)
H- Georeferenciar observaciones meteorológicas
Une 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.
head(observaciones2)## FECHA VALOR STATUS ESTACION year day month
## 7671 01/01/2003 0 M012 2003 1 1
## 7672 02/01/2003 0 M012 2003 2 1
## 7673 03/01/2003 0 M012 2003 3 1
## 7674 04/01/2003 0 M012 2003 4 1
## 7675 05/01/2003 0 M012 2003 5 1
## 7676 06/01/2003 0 M012 2003 6 1head(estaciones2)## CODIGO NOMBRE_DE Lat Long ALTURA PROVINCIA
## 1 M012 LA CUCA -3.492222 -80.05694 53 7
## 2 M032 SANTA ISABEL INAMHI -3.274444 -79.31278 1550 1
## 3 M040 PASAJE -3.329722 -79.78194 40 7
## 4 M142 SARAGURO -3.620556 -79.23222 2525 11
## 5 M179 ARENILLAS -3.560278 -80.05611 60 7
## 6 M180 ZARUMA -3.696944 -79.61611 1100 7observ_ref <- merge(observaciones2, estaciones2,by.x = "ESTACION", by.y = "CODIGO")
class(observ_ref)## [1] "data.frame"head(observ_ref)## ESTACION FECHA VALOR STATUS year day month NOMBRE_DE Lat
## 1 M012 01/01/2003 0 2003 1 1 LA CUCA -3.492222
## 2 M012 02/01/2003 0 2003 2 1 LA CUCA -3.492222
## 3 M012 03/01/2003 0 2003 3 1 LA CUCA -3.492222
## 4 M012 04/01/2003 0 2003 4 1 LA CUCA -3.492222
## 5 M012 05/01/2003 0 2003 5 1 LA CUCA -3.492222
## 6 M012 06/01/2003 0 2003 6 1 LA CUCA -3.492222
## Long ALTURA PROVINCIA coords.x1 coords.x2
## 1 -80.05694 53 7 604749.6 9613947
## 2 -80.05694 53 7 604749.6 9613947
## 3 -80.05694 53 7 604749.6 9613947
## 4 -80.05694 53 7 604749.6 9613947
## 5 -80.05694 53 7 604749.6 9613947
## 6 -80.05694 53 7 604749.6 9613947La variable obtenida es de tipo data.frame, es decir no está georeferenciada. Sin embargo, ahora el data.frame contiene una columna llamada x e y que se utilizarán para georeferenciar con la ayuda de la función coordinates. Luego, se asignará un CRS ya que el SpatialPointsDataFrame no tiene asignado ninguno.
coordinates(observ_ref) <- ~coords.x1+coords.x2 #Opción 1
# coordinates(observ_ref) <- c("coords.x1","coords.x2")#Opción 2
class(observ_ref)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"proj4string(observ_ref) ## [1] NAproj4string(observ_ref) <- CRS(utm17)
proj4string(observ_ref) ## [1] "+proj=utm +zone=17 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0"observ_ref## class : SpatialPointsDataFrame
## features : 6244
## extent : 589069.8, 720854.5, 9561300, 9662797 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=17 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0
## variables : 12
## names : ESTACION, FECHA, VALOR, STATUS, year, day, month, NOMBRE_DE, Lat, Long, ALTURA, PROVINCIA
## min values : M012, 01/01/2003, 0.0, , 2003, 1, 1, CHACRAS, -3.050000, -79.01278, 1100, 1
## max values : M773, 31/12/2003, 91.6, T, 2003, 31, 12, ZARUMA, -3.968056, -80.19806, 60, 7Visualiza los puntos georeferenciados.
l1 = list("sp.points", observ_ref, pch = 19, col="green")
spplot(zona_estudio, "DPA_PROVIN", col="grey", sp.layout = list(l1))
I- Guarda los archivos creados
#Guardar txt
write.table(observaciones2, "observaciones_no_missing.txt")
#Guardar shp
writeOGR(observ_ref, dsn = ".", layer = "observaciones_georef", driver = "ESRI Shapefile",overwrite_layer = TRUE)
#Guardar Geotiff
writeRaster(land_cover2_croped, filename="land_cover2_croped.tif", format="GTiff", overwrite=TRUE)