Raster: álgebra de mapas
Daniela Ballari
Julio, 2018
Maestría en Ecohidrología
Universidad de Cuenca
http://www.ucuenca.edu.ec/maestria-ecohidrologia/
Daniela Ballari (PhD)
dballari@uazuay.edu.ec
Universidad del Azuay
Publicaciones - https://scholar.google.com/citations?user=1VcAm9AAAAAJ
Objetivo de la práctica: Aplicar algebra de mapas local, focal y zonal.
Temario
Parte 1 - Tipos de álgebra de mapas
A-Preparar datos
B-Algebra de mapas local
C-Algebra de mapas focal
D-Algebra de mapas zonal
E-Algebra de mapas global
Parte 2
A-Preparar datos y atributos (valorar criterios)
B-Rasterizar
C-Definir pesos de criterios y aplicar algebra de mapas local con método de ponderación
D-Algebra de mapas zonal: Superponer vías
E-Cuantificar impacto de cada alternativa y representar resultados
¿Cómo hacer esta práctica (Parte 1) con ArcGis? Aqui
¿Cómo hacer esta práctica (Parte 2) con Qgis? Aqui
¿Teoría sobre análisis espacial vectorial?
Olaya Victor (Ed.) (2011) Sistemas de Información Geográfica, Capítulo 14 “Algebra de mapas”
Parte 1 - Tipos de álgebra de mapas
A-Preparar datos
Datos para esta práctica Aqui
library(sf) ## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3library(raster) ## Loading required package: splibrary(latticeExtra)## Loading required package: lattice## Loading required package: RColorBrewersetwd("C:/R_ecohidrologia/Raster") #definir directorio de trabajo Carga las capas de provincias de Ecuador y los 12 rasters de precipitación TRMM para el año 2003
ecuador<-st_read("nxprovincias.shp") #Fuente SNI, archivo shp## Reading layer `nxprovincias' from data source `C:\R_ecohidrologia\Raster\nxprovincias.shp' using driver `ESRI Shapefile'
## Simple feature collection with 25 features and 16 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -732143.5 ymin: 9445216 xmax: 1147852 ymax: 10189400
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defsecuador<-st_transform(ecuador, 4326)
files <- list.files(pattern='.G3.nc') #extraer listado de archivos
files## [1] "TRMM_3B43_ACC.007_20030101_000000.G3.nc"
## [2] "TRMM_3B43_ACC.007_20030201_000000.G3.nc"
## [3] "TRMM_3B43_ACC.007_20030301_000000.G3.nc"
## [4] "TRMM_3B43_ACC.007_20030401_000000.G3.nc"
## [5] "TRMM_3B43_ACC.007_20030501_000000.G3.nc"
## [6] "TRMM_3B43_ACC.007_20030601_000000.G3.nc"
## [7] "TRMM_3B43_ACC.007_20030701_000000.G3.nc"
## [8] "TRMM_3B43_ACC.007_20030801_000000.G3.nc"
## [9] "TRMM_3B43_ACC.007_20030901_000000.G3.nc"
## [10] "TRMM_3B43_ACC.007_20031001_000000.G3.nc"
## [11] "TRMM_3B43_ACC.007_20031101_000000.G3.nc"
## [12] "TRMM_3B43_ACC.007_20031201_000000.G3.nc"Leer archivos
imageraster <- lapply(files,raster)## Loading required namespace: ncdf4spplot(imageraster[[1]], col.regions=colorRampPalette(c("white", "blue"))(255))
Recortar zona de estudio
imageraster <- lapply(imageraster, crop, y=as(ecuador,"Spatial"))#cortar para ecuador
spplot(imageraster[[1]], col.regions=colorRampPalette(c("white", "blue"))(255)) +
as.layer(spplot(as(ecuador[1], "Spatial"), fill="transparent", col="black", under = FALSE))
Establecer máscara de acuerdo a un poligono.
imageraster <- lapply(imageraster, mask, mask=as(ecuador,"Spatial"))#cortar para ecuador
spplot(imageraster[[1]], col.regions=colorRampPalette(c("white", "blue"))(255)) +
as.layer(spplot(as(ecuador[1], "Spatial"), fill="transparent", col="black", under = FALSE))
precipitacion.mensual<- brick(imageraster)
precipitacion.mensual## class : RasterBrick
## dimensions : 27, 67, 1809, 12 (nrow, ncol, ncell, nlayers)
## resolution : 0.25, 0.25 (x, y)
## extent : -92, -75.25, -5, 1.75 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : accumulated_r.1, accumulated_r.2, accumulated_r.3, accumulated_r.4, accumulated_r.5, accumulated_r.6, accumulated_r.7, accumulated_r.8, accumulated_r.9, accumulated_r.10, accumulated_r.11, accumulated_r.12
## min values : 0.03, 1.08, 0.39, 0.00, 0.03, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00
## max values : 439.1327, 495.3195, 415.3514, 526.6675, 551.2563, 523.4909, 386.9567, 320.7042, 382.8305, 404.6749, 427.7702, 500.5126B- Algebra de mapas local
Precipitacion acumulada anual
precipitacion.anual<-sum(precipitacion.mensual)
spplot(precipitacion.anual, col.regions=colorRampPalette(c("white", "blue"))(255)) +
as.layer(spplot(as(ecuador[1], "Spatial"), fill="transparent", col="black", under = FALSE))
precipitacion.anual2<-precipitacion.mensual[[1]]+precipitacion.mensual[[2]]+precipitacion.mensual[[3]]+precipitacion.mensual[[4]]+precipitacion.mensual[[5]]+precipitacion.mensual[[6]]+precipitacion.mensual[[7]]+precipitacion.mensual[[8]]+precipitacion.mensual[[9]]+precipitacion.mensual[[10]]+precipitacion.mensual[[11]]+precipitacion.mensual[[12]]
spplot(precipitacion.anual2, col.regions=colorRampPalette(c("white", "blue"))(255)) +
as.layer(spplot(as(ecuador[1], "Spatial"), fill="transparent", col="black", under = FALSE))
Precipitacion media anual
precipitacion.anual<-mean(precipitacion.mensual)
spplot(precipitacion.anual, col.regions=colorRampPalette(c("white", "blue"))(255)) +
as.layer(spplot(as(ecuador[1], "Spatial"), fill="transparent", col="black", under = FALSE))
C-Algebra de mapas focal
Ventana móvil. Filtro de paso bajo.
precipitacion.filtro <- focal(precipitacion.anual, w=matrix(1,3,3), fun=mean)
spplot(precipitacion.filtro, col.regions=colorRampPalette(c("white", "blue"))(255)) +
as.layer(spplot(as(ecuador[1], "Spatial"), fill="transparent", col="black", under = FALSE))
Ventana móvil. Variabilidad focal
precipitacion.focal.sd <- focal(precipitacion.anual, w=matrix(1,3,3), fun=sd)
spplot(precipitacion.focal.sd, col.regions=colorRampPalette(c("white", "red"))(255)) +
as.layer(spplot(as(ecuador[1], "Spatial"), fill="transparent", col="black", under = FALSE))
D- Algebra de mapas zonal
zonal <- extract(precipitacion.anual, as(ecuador, "Spatial"), fun=mean, na.rm=TRUE, df=TRUE)
zonal## ID layer
## 1 1 76.475982
## 2 2 107.345040
## 3 3 78.186558
## 4 4 172.317750
## 5 5 109.514125
## 6 6 99.993900
## 7 7 38.606768
## 8 8 171.639052
## 9 9 63.907826
## 10 10 114.935725
## 11 11 53.845399
## 12 12 131.689340
## 13 13 87.534196
## 14 14 191.998192
## 15 15 233.148746
## 16 16 299.467908
## 17 17 112.898620
## 18 18 125.161280
## 19 19 116.863880
## 20 20 1.558851
## 21 21 246.191826
## 22 22 283.112023
## 23 23 169.163751
## 24 24 38.184608
## 25 25 172.383794E- Algebra de mapas global
Precipitacion anual total en el Ecuador
precipitacion.mensual.global<- cellStats(precipitacion.mensual,sum)
barplot(precipitacion.mensual.global)
precipitation.total<- cellStats(precipitacion.anual,sum)
precipitation.total## [1] 54178.52precipitation.media<- cellStats(precipitacion.anual,mean)
precipitation.media## [1] 163.6813Parte 2 - Ejemplo de álgebra de mapas para toma de decisiones espaciales
A-Preparar datos y atributos (valorar criterios)
Cargar datos
prov_Guayas<-st_read("Prov_Guayas.shp") #Fuente SNI, archivo shp## Reading layer `Prov_Guayas' from data source `C:\R_ecohidrologia\Raster\Prov_Guayas.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 548698.8 ymin: 9661256 xmax: 711084.4 ymax: 9907510
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defsvias<-st_read("vias.shp") #Fuente elaboración propia, archivo shp## Reading layer `vias' from data source `C:\R_ecohidrologia\Raster\vias.shp' using driver `ESRI Shapefile'
## Simple feature collection with 3 features and 1 field
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 563604.6 ymin: 9678370 xmax: 687246.3 ymax: 9895096
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defssummary(vias)#explorar contenido de datos## id geometry
## Min. :1.0 LINESTRING :3
## 1st Qu.:1.5 epsg:32717 :0
## Median :2.0 +proj=utm ...:0
## Mean :2.0
## 3rd Qu.:2.5
## Max. :3.0plot(st_geometry(prov_Guayas),col="grey")
plot(st_geometry(vias),col="red", add=TRUE)
clima_Guayas<-st_read("Clima_ Guayas.shp") #Fuente SNI, archivo shp## Reading layer `Clima_ Guayas' from data source `C:\R_ecohidrologia\Raster\Clima_ Guayas.shp' using driver `ESRI Shapefile'
## Simple feature collection with 64 features and 10 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 548698.8 ymin: 9661256 xmax: 711084.4 ymax: 9907510
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defssummary(clima_Guayas)#explorar contenido de datos## codigo descripcio area_ha
## Ab:44 ECUATORIAL MESOTERMICO SEMI-HUMEDO: 1 Min. : 12
## Ah: 6 TROPICAL MEGATERMICO HUMEDO : 1 1st Qu.: 140
## Ar:12 TROPICAL MEGATERMICO SECO :44 Median : 515
## Aw: 1 TROPICAL MEGATERMICO SEMI ARIDO :12 Mean : 172170
## Ch: 1 TROPICAL MEGATERMICO SEMI HUMEDO : 6 3rd Qu.: 2331
## Max. :3873572
## DPA_PROVIN DPA_DESPRO DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO
## 09:64 GUAYAS:64 Min. :0 2011:64 03:64 02:64
## 1st Qu.:0
## Median :0
## Mean :0
## 3rd Qu.:0
## Max. :0
## PEE_CODIGO geometry
## 593:64 MULTIPOLYGON :64
## epsg:32717 : 0
## +proj=utm ...: 0
##
##
## plot(clima_Guayas["descripcio"])
hidrogeologia_Guayas<-st_read("Hidrogeologia_Guayas.shp") #Fuente SNI, archivo shp## Reading layer `Hidrogeologia_Guayas' from data source `C:\R_ecohidrologia\Raster\Hidrogeologia_Guayas.shp' using driver `ESRI Shapefile'
## Simple feature collection with 760 features and 16 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 548698.8 ymin: 9661256 xmax: 711084.4 ymax: 9907510
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defssummary(hidrogeologia_Guayas)#explorar contenido de datos## ID SIMBOLO PERM
## Min. : 0 SALITRALES:238 A4 :338
## 1st Qu.: 2450 Da : 65 Wn :130
## Median : 5306 Lago : 47 B3 : 97
## Mean : 6730 Daem : 42 A3 : 88
## 3rd Qu.:11104 Mar : 37 A1 : 65
## Max. :11422 Pl DB : 36 A2 : 21
## (Other) :295 (Other): 21
## LITOLOGIA FORMACION1
## SALITRALES :182 DEPOSITO ALUVIAL DE ESTERO: 25
## DEPOSITO ALUVIAL DE ESTERO :103 F.BORBON : 23
## DEPOSITO ALUVIAL : 31 F.CAYO 2500 M : 21
## DEPOSITO ALUVIAL: ARCILLAS, ARENAS: 28 F.ONZOLE : 19
## ARCILLAS, LIMOS, ARENISCAS : 25 U.MACUCHI : 16
## (Other) :261 (Other) :120
## NA's :130 NA's :536
## EDAD PERMIABILI
## Cuaternaria :431 BAJA :338
## Cretaceo : 48 MUY BAJA : 97
## Plioceno : 38 MEDIA : 88
## Mioceno : 23 GENERALMENTE ALTA: 65
## Mio-Plioceno: 19 MEDIA A ALTA : 21
## (Other) : 58 (Other) : 21
## NA's :143 NA's :130
## TIPO_PERM
## FISURACION :101
## POROSIDAD INTERGRANULAR :512
## POROSIDAD INTERGRANULAR Y FISURACION-ROCAS SIN IMPORTANCIA HIDROGEOLOGICA: 17
## NA's :130
##
##
##
## SIMBOLO1 DPA_PROVIN DPA_DESPRO DPA_VALOR DPA_ANIO
## SALITRALES:238 09:760 GUAYAS:760 Min. :0 2011:760
## Wn : 76 1st Qu.:0
## Da : 65 Median :0
## Daem : 42 Mean :0
## Mar : 37 3rd Qu.:0
## Pl DB : 36 Max. :0
## (Other) :266
## REI_CODIGO REN_CODIGO PEE_CODIGO geometry
## 03:760 02:760 593:760 MULTIPOLYGON :760
## epsg:32717 : 0
## +proj=utm ...: 0
##
##
##
## plot(hidrogeologia_Guayas["PERMIABILI"])
zonas_inundables_Guayas<-st_read("Zonas_Inundables_Guayas.shp") #Fuente SNI, archivo shp## Reading layer `Zonas_Inundables_Guayas' from data source `C:\R_ecohidrologia\Raster\Zonas_Inundables_Guayas.shp' using driver `ESRI Shapefile'
## Simple feature collection with 37 features and 13 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 611744.9 ymin: 9680116 xmax: 651935.2 ymax: 9867700
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defssummary(zonas_inundables_Guayas)#explorar contenido de datos## fcode
## BH092:37
##
##
##
##
##
## descripcio
## TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSI<d3>N DE LAS AGUAS DE MAREA:37
##
##
##
##
##
## txt SHAPE_Leng Shape_Le_1 Shape_Area
## T. S. I.:33 Min. : 305.9 Min. : 305.9 Min. : 4211
## NA's : 4 1st Qu.: 1185.1 1st Qu.: 1185.1 1st Qu.: 44077
## Median : 2240.6 Median : 2240.6 Median : 94380
## Mean : 3711.8 Mean : 3711.8 Mean : 932462
## 3rd Qu.: 4513.9 3rd Qu.: 4513.9 3rd Qu.: 160755
## Max. :22837.7 Max. :22837.7 Max. :26066129
## DPA_PROVIN DPA_DESPRO DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO
## 09:37 GUAYAS:37 Min. :0 2011:37 03:37 02:37
## 1st Qu.:0
## Median :0
## Mean :0
## 3rd Qu.:0
## Max. :0
## PEE_CODIGO geometry
## 593:37 POLYGON :37
## epsg:32717 : 0
## +proj=utm ...: 0
##
##
## plot(st_geometry(prov_Guayas),col="grey")
plot(st_geometry(zonas_inundables_Guayas),col="blue", border="blue", add=TRUE)
Valorar atributos
#CLIMA
names(clima_Guayas)# explorar nombres de campos## [1] "codigo" "descripcio" "area_ha" "DPA_PROVIN" "DPA_DESPRO"
## [6] "DPA_VALOR" "DPA_ANIO" "REI_CODIGO" "REN_CODIGO" "PEE_CODIGO"
## [11] "geometry"head(clima_Guayas)# explorar unos pocos elementos de la tabla asociada## Simple feature collection with 6 features and 10 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 548698.8 ymin: 9714390 xmax: 710783.2 ymax: 9824784
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## codigo descripcio area_ha DPA_PROVIN
## 1 Aw TROPICAL MEGATERMICO HUMEDO 1.981810e+06 09
## 2 Ab TROPICAL MEGATERMICO SECO 9.899287e+05 09
## 3 Ah TROPICAL MEGATERMICO SEMI HUMEDO 2.318054e+03 09
## 4 Ah TROPICAL MEGATERMICO SEMI HUMEDO 2.185124e+03 09
## 5 Ah TROPICAL MEGATERMICO SEMI HUMEDO 5.957881e+01 09
## 6 Ah TROPICAL MEGATERMICO SEMI HUMEDO 1.707123e+02 09
## DPA_DESPRO DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO PEE_CODIGO
## 1 GUAYAS 0 2011 03 02 593
## 2 GUAYAS 0 2011 03 02 593
## 3 GUAYAS 0 2011 03 02 593
## 4 GUAYAS 0 2011 03 02 593
## 5 GUAYAS 0 2011 03 02 593
## 6 GUAYAS 0 2011 03 02 593
## geometry
## 1 MULTIPOLYGON (((710570.1 97...
## 2 MULTIPOLYGON (((565511.5 97...
## 3 MULTIPOLYGON (((626057.6 97...
## 4 MULTIPOLYGON (((629105.9 97...
## 5 MULTIPOLYGON (((630022.2 97...
## 6 MULTIPOLYGON (((619764.7 97...table(clima_Guayas$descripcio) # resumen de valores y frecuencias del 'campo descripcio'##
## ECUATORIAL MESOTERMICO SEMI-HUMEDO TROPICAL MEGATERMICO HUMEDO
## 1 1
## TROPICAL MEGATERMICO SECO TROPICAL MEGATERMICO SEMI ARIDO
## 44 12
## TROPICAL MEGATERMICO SEMI HUMEDO
## 6# 1- Definir escala de valoración: en este caso será de 1 a 5, con 1
#los que produzca mayor impacto para la via y 5 lo que produzca menor impacto.
# 2- Asignar a cada atributo la valoración. Este paso es debido
#a que se necesita contar con valores numéricos para la rasterizacion.
#Luego cada valor será una clase en el archivo raster.
clima_Guayas$reclass[clima_Guayas$descripcio =="TROPICAL MEGATERMICO SECO"] = 5
clima_Guayas$reclass[clima_Guayas$descripcio =="TROPICAL MEGATERMICO SEMI ARIDO"] = 4
clima_Guayas$reclass[clima_Guayas$descripcio =="ECUATORIAL MESOTERMICO SEMI-HUMEDO"] = 2
clima_Guayas$reclass[clima_Guayas$descripcio =="TROPICAL MEGATERMICO SEMI HUMEDO"] = 2
clima_Guayas$reclass[clima_Guayas$descripcio =="TROPICAL MEGATERMICO HUMEDO"] = 1
table(clima_Guayas$reclass)# nuevos valores y sus frecuencias##
## 1 2 4 5
## 1 7 12 44#HIDROLOGIA
names(hidrogeologia_Guayas)# explorar nombres de campos## [1] "ID" "SIMBOLO" "PERM" "LITOLOGIA" "FORMACION1"
## [6] "EDAD" "PERMIABILI" "TIPO_PERM" "SIMBOLO1" "DPA_PROVIN"
## [11] "DPA_DESPRO" "DPA_VALOR" "DPA_ANIO" "REI_CODIGO" "REN_CODIGO"
## [16] "PEE_CODIGO" "geometry"head(hidrogeologia_Guayas) # explorar unos pocos elementos de la tabla asociada## Simple feature collection with 6 features and 16 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 580208.5 ymin: 9821276 xmax: 609708.2 ymax: 9854039
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## ID SIMBOLO PERM LITOLOGIA FORMACION1
## 1 6 Pl DB A3 ARENISCA DELESNABLE CON BANCOS CALCAREOS F.BORBON
## 2 9 Da A1 DEPOSITO ALUVIAL <NA>
## 3 26 MP l DO B3 ARCILLA, LUTITA F.ONZOLE
## 4 32 Pl DB A3 ARENISCA DELESNABLE CON BANCOS CALCAREOS F.BORBON
## 5 34 Pl DB A3 ARENISCA DELESNABLE CON BANCOS CALCAREOS F.BORBON
## 6 37 Lago Wn <NA> <NA>
## EDAD PERMIABILI TIPO_PERM SIMBOLO1
## 1 Plioceno MEDIA POROSIDAD INTERGRANULAR Pl DB
## 2 Cuaternaria GENERALMENTE ALTA POROSIDAD INTERGRANULAR Da
## 3 Mio-Plioceno MUY BAJA FISURACION MP l DO
## 4 Plioceno MEDIA POROSIDAD INTERGRANULAR Pl DB
## 5 Plioceno MEDIA POROSIDAD INTERGRANULAR Pl DB
## 6 <NA> <NA> <NA> Wn
## DPA_PROVIN DPA_DESPRO DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO
## 1 09 GUAYAS 0 2011 03 02
## 2 09 GUAYAS 0 2011 03 02
## 3 09 GUAYAS 0 2011 03 02
## 4 09 GUAYAS 0 2011 03 02
## 5 09 GUAYAS 0 2011 03 02
## 6 09 GUAYAS 0 2011 03 02
## PEE_CODIGO geometry
## 1 593 MULTIPOLYGON (((592144.7 98...
## 2 593 MULTIPOLYGON (((605767.8 98...
## 3 593 MULTIPOLYGON (((582842.6 98...
## 4 593 MULTIPOLYGON (((588295 9834...
## 5 593 MULTIPOLYGON (((599285.9 98...
## 6 593 MULTIPOLYGON (((598830.9 98...table(hidrogeologia_Guayas$PERMIABILI) # resumen de valores y frecuencias del 'campo descripcio'##
## BAJA GENERALMENTE ALTA
## 338 65
## GENERALMENTE BAJA MEDIA
## 4 88
## MEDIA A ALTA MUY BAJA
## 21 97
## PRACTICAMENTE IMPERMEABLE
## 17hidrogeologia_Guayas$reclass[hidrogeologia_Guayas$PERMIABILI =="PRACTICAMENTE IMPERMEABLE"] = 5
hidrogeologia_Guayas$reclass[hidrogeologia_Guayas$PERMIABILI =="MUY BAJA"] = 4
hidrogeologia_Guayas$reclass[hidrogeologia_Guayas$PERMIABILI =="BAJA"] = 3
hidrogeologia_Guayas$reclass[hidrogeologia_Guayas$PERMIABILI =="GENERALMENTE BAJA"] = 2
hidrogeologia_Guayas$reclass[hidrogeologia_Guayas$PERMIABILI =="MEDIA"] = 2
hidrogeologia_Guayas$reclass[hidrogeologia_Guayas$PERMIABILI =="MEDIA A ALTA"] = 1
hidrogeologia_Guayas$reclass[hidrogeologia_Guayas$PERMIABILI =="GENERALMENTE ALTA"] = 1
table(hidrogeologia_Guayas$reclass)##
## 1 2 3 4 5
## 86 92 338 97 17#ZONAS INUNDABLES
names(zonas_inundables_Guayas)# explorar nombres de campos## [1] "fcode" "descripcio" "txt" "SHAPE_Leng" "Shape_Le_1"
## [6] "Shape_Area" "DPA_PROVIN" "DPA_DESPRO" "DPA_VALOR" "DPA_ANIO"
## [11] "REI_CODIGO" "REN_CODIGO" "PEE_CODIGO" "geometry"head(zonas_inundables_Guayas) # explorar unos pocos elementos de la tabla asociada## Simple feature collection with 6 features and 13 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 613906.1 ymin: 9846156 xmax: 623700.1 ymax: 9856021
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## fcode
## 1 BH092
## 2 BH092
## 3 BH092
## 4 BH092
## 5 BH092
## 6 BH092
## descripcio
## 1 TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSI<d3>N DE LAS AGUAS DE MAREA
## 2 TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSI<d3>N DE LAS AGUAS DE MAREA
## 3 TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSI<d3>N DE LAS AGUAS DE MAREA
## 4 TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSI<d3>N DE LAS AGUAS DE MAREA
## 5 TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSI<d3>N DE LAS AGUAS DE MAREA
## 6 TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSI<d3>N DE LAS AGUAS DE MAREA
## txt SHAPE_Leng Shape_Le_1 Shape_Area DPA_PROVIN DPA_DESPRO
## 1 T. S. I. 2154.2896 2154.2896 77375.229 09 GUAYAS
## 2 T. S. I. 639.9910 639.9910 23492.236 09 GUAYAS
## 3 T. S. I. 416.3941 416.3941 7622.892 09 GUAYAS
## 4 T. S. I. 2240.5962 2240.5962 55050.758 09 GUAYAS
## 5 <NA> 2010.2875 2010.2875 141961.969 09 GUAYAS
## 6 <NA> 2427.1830 2427.1830 130237.668 09 GUAYAS
## DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO PEE_CODIGO
## 1 0 2011 03 02 593
## 2 0 2011 03 02 593
## 3 0 2011 03 02 593
## 4 0 2011 03 02 593
## 5 0 2011 03 02 593
## 6 0 2011 03 02 593
## geometry
## 1 POLYGON ((617613.7 9853969,...
## 2 POLYGON ((623654.9 9855808,...
## 3 POLYGON ((614740.1 9854033,...
## 4 POLYGON ((614215 9853555, 6...
## 5 POLYGON ((619999.5 9846900,...
## 6 POLYGON ((620162.7 9846649,...table(zonas_inundables_Guayas$descripcio) # resumen de valores y frecuencias del 'campo descripcio'##
## TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSI<d3>N DE LAS AGUAS DE MAREA
## 37zonas_inundables_Guayas$reclass = 1
table(zonas_inundables_Guayas$reclass)##
## 1
## 37B-Rasterizar
r <- raster(as(prov_Guayas, "Spatial"), ncols = 1000, nrows = 1000)
#definir raster (vacío)
extent(r) ## class : Extent
## xmin : 548698.8
## xmax : 711084.4
## ymin : 9661256
## ymax : 9907510extent(prov_Guayas) #toma la extensión espacial de capa existente## class : Extent
## xmin : 548698.8
## xmax : 711084.4
## ymin : 9661256
## ymax : 9907510projection(r) ## [1] "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"st_crs(prov_Guayas) #toma CRS de capa existente## Coordinate Reference System:
## EPSG: 32717
## proj4string: "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs"res(r)<-c(500,500) #asignar resolucion espacialclima_ras<- rasterize(clima_Guayas, r, "reclass", fun='min')
# rasterize(archivo vectorial, raster base, "campo numerico para rasterizar", fun='funcion a aplicar por la rasterización')
hidrogeologia_ras<- rasterize(hidrogeologia_Guayas, r, "reclass", fun='min')
zonas_inundables_ras<- rasterize(zonas_inundables_Guayas, r, "reclass", fun='min')#Durante la rasterización las zonas fuera de los poligons son asingados con NA.
head(is.na(zonas_inundables_ras))## 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 1 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 2 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 3 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 4 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 5 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 6 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 7 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 8 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 9 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 10 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 15 16 17 18 19 20
## 1 TRUE TRUE TRUE TRUE TRUE TRUE
## 2 TRUE TRUE TRUE TRUE TRUE TRUE
## 3 TRUE TRUE TRUE TRUE TRUE TRUE
## 4 TRUE TRUE TRUE TRUE TRUE TRUE
## 5 TRUE TRUE TRUE TRUE TRUE TRUE
## 6 TRUE TRUE TRUE TRUE TRUE TRUE
## 7 TRUE TRUE TRUE TRUE TRUE TRUE
## 8 TRUE TRUE TRUE TRUE TRUE TRUE
## 9 TRUE TRUE TRUE TRUE TRUE TRUE
## 10 TRUE TRUE TRUE TRUE TRUE TRUE#Aqui se convierten los NA en 0 para posibitar el algebra de mapas posterior.
zonas_inundables_ras[is.na(zonas_inundables_ras)] <- 0 Visualizar
plot(clima_Guayas["descripcio"],main="clima vector", col=rev(heat.colors(5)))## Warning in plot.sf(clima_Guayas["descripcio"], main = "clima vector", col
## = rev(heat.colors(5))): col is not of length 1 or ncol(x): colors will be
## recycled; use pal to specify a color palette
plot(clima_ras, main="clima raster", col=rev(heat.colors(20)))
plot(hidrogeologia_Guayas["PERMIABILI"], main="Hidrogeología vector",
col=rev(heat.colors(7)))## Warning in plot.sf(hidrogeologia_Guayas["PERMIABILI"], main =
## "Hidrogeología vector", : col is not of length 1 or ncol(x): colors will be
## recycled; use pal to specify a color palette
plot(hidrogeologia_ras, main="hidrogeologia raster",
col=rev(heat.colors(20)))
plot(st_geometry(prov_Guayas),border="white", col="grey")
plot(st_geometry(zonas_inundables_Guayas),col="blue", border="blue", add=TRUE)
plot(zonas_inundables_ras, main="area inundacion raster",
col=c("grey", "blue"))
Explorar funciones de libreria raster
#values(clima_ras)
projection(clima_ras)## [1] "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"ncell(clima_ras)## [1] 160225clima_ras[250, 250]##
## 2cellStats(clima_ras,max) ## [1] 5cellStats(clima_ras,min) ## [1] 1cellStats(clima_ras,mean) ## [1] 3.268329cellStats(clima_ras,sd) ## [1] 1.452143C-Definir pesos de criterios y aplicar algebra de mapas local mediante el método de ponderación
Definir pesos de criterios + clima=0.2 + hidrogeologia=0.5 + inundacion=0.3
algebra_mapas <- (0.2*clima_ras) + (0.5*hidrogeologia_ras) + (0.3*zonas_inundables_ras)
algebra_mapas## class : RasterLayer
## dimensions : 493, 325, 160225 (nrow, ncol, ncell)
## resolution : 500, 500 (x, y)
## extent : 548698.8, 711198.8, 9661010, 9907510 (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 : in memory
## names : layer
## values : 0.7, 3.5 (min, max)plot(algebra_mapas, col=rev(heat.colors(20)))
cellStats(algebra_mapas,max) ## [1] 3.5cellStats(algebra_mapas,min) ## [1] 0.7cellStats(algebra_mapas,mean) ## [1] 1.852662D-Algebra de mapas zonal: Superponer vías
extract <- extract(algebra_mapas,vias)
class(extract)## [1] "list"summary(extract)## Length Class Mode
## [1,] 481 -none- numeric
## [2,] 647 -none- numeric
## [3,] 532 -none- numericE-Cuantificar impacto de cada alternativa y representar resultados
Alternativa 1
summary(extract[[1]])## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.900 0.900 1.400 1.387 1.900 1.900 4length(extract[[1]])## [1] 481sum(extract[[1]], na.rm = TRUE)## [1] 661.8mean(extract[[1]], na.rm = TRUE)## [1] 1.387421sd(extract[[1]], na.rm = TRUE)## [1] 0.4650032Alternativa 2
summary(extract[[2]])## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.900 0.900 1.400 1.423 2.000 3.000 9length(extract[[2]])## [1] 647sum(extract[[2]], na.rm = TRUE)## [1] 908mean(extract[[2]], na.rm = TRUE)## [1] 1.423197sd(extract[[2]], na.rm = TRUE)## [1] 0.6014852Alternativa 3
summary(extract[[3]])## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.900 0.900 1.900 1.651 1.900 2.500 1length(extract[[3]])## [1] 532sum(extract[[3]], na.rm = TRUE)## [1] 876.6mean(extract[[3]], na.rm = TRUE)## [1] 1.650847sd(extract[[3]], na.rm = TRUE)## [1] 0.5874266sum <- c(sum(extract[[1]], na.rm = TRUE), sum(extract[[2]],
na.rm = TRUE), sum(extract[[3]], na.rm = TRUE))two.par <- par(mfrow=c(1, 2))
par(mfrow=c(1, 2))
boxplot(extract, main="Boxplot")
barplot(sum, main="valoracion total (suma)", xlab="alternativa")
¿Qué alternativa de vía eligirías?
Trabajo Calificado 6p
- Calcula la precipitación promedio para todo el Ecuador (Álgebra de mapas local). Reporta una captura de pantalla con simbología adecuada para la precipitación promedio. 0.2p
- Calcula la precipitación promedio para una provincia del Ecuador (Algebra de mapas zonal). Reporta una captura de pantalla con la precipitación promedio para tu provincia. 0.2p
- Recorta la capa resultante con tu provincia. Reporta una captura de pantalla del mapa recortado de la precipitación promedio en la provincia asignada. 0.2p
Has llegado al final del material!