Lectura de archivos de geoinformación y sistemas de referencia espacial
Daniela Ballari
Octubre, 2015
Maestría en Ecohidrologí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
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 metereológicas
I- Guarda los archivos creados
A- Librerías
Para comenzar carga las gdal,sp, rastery rgeos y define el directorio de trabajo.
library(sp)
library(rgdal)
library(raster)
library(rgeos)
setwd("C:/R_ecohidrologia/SRS")Explora algunas de las opciones de la libreria rgdal.
ogrDrivers() # Listado de 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 FALSElist.files() # Listado de archivos en el directorio de trabajo## [1] "elevacion.tif" "elevacion.tif.aux.xml"
## [3] "estaciones_32717.dbf" "estaciones_32717.prj"
## [5] "estaciones_32717.qpj" "estaciones_32717.shp"
## [7] "estaciones_32717.shx" "land_cover2_croped.tif"
## [9] "LC_5min_global_2012.tif" "lectura_srs2.html"
## [11] "lectura_srs2.Rmd" "observaciones_2003.csv"
## [13] "observaciones_georef.dbf" "observaciones_georef.prj"
## [15] "observaciones_georef.shp" "observaciones_georef.shx"
## [17] "observaciones_no_missing.txt" "prov_zona_estudio.dbf"
## [19] "prov_zona_estudio.prj" "prov_zona_estudio.shp"
## [21] "prov_zona_estudio.shx" "rio_32717_zona_estudio.dbf"
## [23] "rio_32717_zona_estudio.prj" "rio_32717_zona_estudio.qpj"
## [25] "rio_32717_zona_estudio.shp" "rio_32717_zona_estudio.shx"
## [27] "rio_torrente.dbf" "rio_torrente.prj"
## [29] "rio_torrente.shp" "rio_torrente.shx"ogrListLayers(".") #Archivos georreferenciados shp en el directorio de trabajo## [1] "estaciones_32717" "observaciones_georef"
## [3] "prov_zona_estudio" "rio_32717_zona_estudio"
## [5] "rio_torrente"
## attr(,"driver")
## [1] "ESRI Shapefile"
## attr(,"nlayers")
## [1] 5ogrInfo(".","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 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 IntegerB- Cargar 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 de datos: editado INEC## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "prov_zona_estudio"
## with 7 features
## It has 11 fieldsC- 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. En cambio, si se quisiera cambiar el sistema de referencia (transormar) 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 1 #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 spTransform in package rgdalzona_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, 19Visualiza 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Ã<U+0091>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"D- Carga una capa shp de puntos
estaciones<-readOGR(".","estaciones_32717") #estaciones metereológicas## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "estaciones_32717"
## with 40 features
## It has 14 fieldsproj4string(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 Min. :170.0
## 2P - 3P: 2 M032 : 1 BALSAS : 1 CO: 9 1st Qu.:180.0
## 3MP :15 M040 : 1 CHACRAS : 1 CP: 2 Median :195.0
## 3P :21 M142 : 1 COCHAPAMBA-QUINGEO: 1 PG: 2 Mean :205.2
## M179 : 1 CUMBE : 1 PV:25 3rd Qu.:200.0
## M180 : 1 EL SALADO-PREDESUR: 1 Max. :280.0
## (Other):34 (Other) :34
## Lat Long ALTURA PROVINCIA
## Min. :-3.968 Min. :-80.22 Min. : 5.0 Min. : 1.0
## 1st Qu.:-3.724 1st Qu.:-79.77 1st Qu.: 301.2 1st Qu.: 1.0
## Median :-3.478 Median :-79.62 Median :1113.0 Median : 7.0
## Mean :-3.492 Mean :-79.55 Mean :1271.9 Mean : 5.9
## 3rd Qu.:-3.284 3rd Qu.:-79.25 3rd Qu.:2335.0 3rd Qu.: 7.0
## Max. :-3.004 Max. :-78.92 Max. :3450.0 Max. :19.0
##
## INSTITUCIO FECHA_DE_I FECHA_DE_L x
## INAMHI :23 1972 / 1 / 1 : 3 1958 / 4 / 1 : 1 Min. :586178
## INECEL : 4 1981 / 2 / 1 : 3 1966 / 2 / 2 : 1 1st Qu.:636529
## INERHI : 1 1963 / 5 / 16: 2 1974 / 1 / 15: 1 Median :653549
## PREDESUR:12 1973 / 2 / 1 : 2 1984 / 1 / 1 : 1 Mean :661292
## 1978 / 8 / 7 : 2 NA's :36 3rd Qu.:694320
## (Other) :25 Max. :731094
## NA's : 3
## y
## Min. :9561300
## 1st Qu.:9588266
## Median :9615409
## Mean :9613894
## 3rd Qu.:9636876
## Max. :9667788
## Selecciona determinadas columnas del SpatialPointsDataFrame.
estaciones2<- estaciones[,c(2,3,6,7,8,9)]#opción 1: 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
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.
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.
####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)
#spplot(land_cover) # funciona, pero muy lento
####Opción 1 raster
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 : C:\R_ecohidrologia\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- Georreferenciar observaciones metereoló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á 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) <- ~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, 5, 1
## max values : M773, 31/12/2003, 91.6, T, 2003, 31, 12, ZARUMA, -3.968056, -80.19806, 2750, 19Visualiza 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)Has llegado al final de este material!