Raster: álgebra de mapas
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
Objetivo de la práctica: Aplicar algebra de mapas y toma de decisiones para seleccionar la mejor de tres alternativas de vias de acuerdo a factores ambientales. Se utilizarán tres capas de información: clima, hidrogeología y areas inundables.
Temario
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 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” — ###A-Preparar datos y atributos (valorar criterios)
library(rgdal)
library(raster)
setwd("C:/R_ecohidrologia/Raster") #definir directorio de trabajo #cargar datos
prov_Guayas<-readOGR(".","Prov_Guayas") #Fuente SNI, archivo shp## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "Prov_Guayas"
## with 1 features
## It has 7 fieldsvias<-readOGR(".","vias") #Fuente elaboración propia, archivo shp## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "vias"
## with 3 features
## It has 1 fieldssummary(vias)#explorar contenido de datos## Object of class SpatialLinesDataFrame
## Coordinates:
## min max
## x 563604.6 687246.3
## y 9678370.5 9895096.1
## Is projected: TRUE
## proj4string :
## [+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## id
## Min. :1.0
## 1st Qu.:1.5
## Median :2.0
## Mean :2.0
## 3rd Qu.:2.5
## Max. :3.0lines = list("sp.lines", vias, col="red")
spplot(prov_Guayas, "DPA_DESPRO",col="white", fill="grey", sp.layout = list(lines),colorkey=FALSE)
clima_Guayas<-readOGR(".","Clima_ Guayas") #Fuente SNI, archivo shp## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "Clima_ Guayas"
## with 64 features
## It has 10 fieldssummary(clima_Guayas)#explorar contenido de datos## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 548698.8 711084.4
## y 9661255.8 9907509.5
## Is projected: TRUE
## proj4string :
## [+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## 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
## 593:64
##
##
##
##
## spplot(clima_Guayas, "descripcio")
hidrogeologia_Guayas<-readOGR(".","Hidrogeologia_Guayas") #Fuente SNI, archivo shp## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "Hidrogeologia_Guayas"
## with 760 features
## It has 16 fieldssummary(hidrogeologia_Guayas)#explorar contenido de datos## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 548698.8 711084.4
## y 9661255.8 9907509.5
## Is projected: TRUE
## proj4string :
## [+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## 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
## 03:760 02:760 593:760
##
##
##
##
##
## spplot(hidrogeologia_Guayas, "PERMIABILI")
zonas_inundables_Guayas<-readOGR(".","Zonas_Inundables_Guayas") #Fuente SNI, archivo shp## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "Zonas_Inundables_Guayas"
## with 37 features
## It has 13 fieldssummary(zonas_inundables_Guayas)#explorar contenido de datos## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 611744.9 651935.2
## y 9680115.8 9867699.9
## Is projected: TRUE
## proj4string :
## [+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## fcode
## BH092:37
##
##
##
##
##
## descripcio
## TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSIÓ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
## 593:37
##
##
##
##
## polygon = list("sp.polygons", zonas_inundables_Guayas, col="blue", fill="blue")
spplot(prov_Guayas, "DPA_DESPRO",col="white", fill="grey", sp.layout = list(polygon),colorkey=FALSE)
#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"head(clima_Guayas)# explorar unos pocos elementos de la tabla asociada## codigo descripcio area_ha DPA_PROVIN
## 0 Aw TROPICAL MEGATERMICO HUMEDO 1.981810e+06 09
## 1 Ab TROPICAL MEGATERMICO SECO 9.899287e+05 09
## 2 Ah TROPICAL MEGATERMICO SEMI HUMEDO 2.318054e+03 09
## 3 Ah TROPICAL MEGATERMICO SEMI HUMEDO 2.185124e+03 09
## 4 Ah TROPICAL MEGATERMICO SEMI HUMEDO 5.957881e+01 09
## 5 Ah TROPICAL MEGATERMICO SEMI HUMEDO 1.707123e+02 09
## DPA_DESPRO DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO PEE_CODIGO
## 0 GUAYAS 0 2011 03 02 593
## 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 593table(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"head(hidrogeologia_Guayas) # explorar unos pocos elementos de la tabla asociada## ID SIMBOLO PERM LITOLOGIA FORMACION1
## 0 6 Pl DB A3 ARENISCA DELESNABLE CON BANCOS CALCAREOS F.BORBON
## 1 9 Da A1 DEPOSITO ALUVIAL <NA>
## 2 26 MP l DO B3 ARCILLA, LUTITA F.ONZOLE
## 3 32 Pl DB A3 ARENISCA DELESNABLE CON BANCOS CALCAREOS F.BORBON
## 4 34 Pl DB A3 ARENISCA DELESNABLE CON BANCOS CALCAREOS F.BORBON
## 5 37 Lago Wn <NA> <NA>
## EDAD PERMIABILI TIPO_PERM SIMBOLO1
## 0 Plioceno MEDIA POROSIDAD INTERGRANULAR Pl DB
## 1 Cuaternaria GENERALMENTE ALTA POROSIDAD INTERGRANULAR Da
## 2 Mio-Plioceno MUY BAJA FISURACION MP l DO
## 3 Plioceno MEDIA POROSIDAD INTERGRANULAR Pl DB
## 4 Plioceno MEDIA POROSIDAD INTERGRANULAR Pl DB
## 5 <NA> <NA> <NA> Wn
## DPA_PROVIN DPA_DESPRO DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO
## 0 09 GUAYAS 0 2011 03 02
## 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
## PEE_CODIGO
## 0 593
## 1 593
## 2 593
## 3 593
## 4 593
## 5 593table(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"head(zonas_inundables_Guayas) # explorar unos pocos elementos de la tabla asociada## fcode
## 0 BH092
## 1 BH092
## 2 BH092
## 3 BH092
## 4 BH092
## 5 BH092
## descripcio
## 0 TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSIÓN DE LAS AGUAS DE MAREA
## 1 TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSIÓN DE LAS AGUAS DE MAREA
## 2 TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSIÓN DE LAS AGUAS DE MAREA
## 3 TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSIÓN DE LAS AGUAS DE MAREA
## 4 TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSIÓN DE LAS AGUAS DE MAREA
## 5 TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSIÓN DE LAS AGUAS DE MAREA
## txt SHAPE_Leng Shape_Le_1 Shape_Area DPA_PROVIN DPA_DESPRO
## 0 T. S. I. 2154.2896 2154.2896 77375.229 09 GUAYAS
## 1 T. S. I. 639.9910 639.9910 23492.236 09 GUAYAS
## 2 T. S. I. 416.3941 416.3941 7622.892 09 GUAYAS
## 3 T. S. I. 2240.5962 2240.5962 55050.758 09 GUAYAS
## 4 <NA> 2010.2875 2010.2875 141961.969 09 GUAYAS
## 5 <NA> 2427.1830 2427.1830 130237.668 09 GUAYAS
## DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO PEE_CODIGO
## 0 0 2011 03 02 593
## 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 593table(zonas_inundables_Guayas$descripcio) # resumen de valores y frecuencias del 'campo descripcio'##
## TRAMO QUE ESTA EXCEPCIONALMENTE CUBIERTO POR AGUA, CON EXCLUSIÓN DE LAS AGUAS DE MAREA
## 37zonas_inundables_Guayas$reclass = 1
table(zonas_inundables_Guayas$reclass)##
## 1
## 37B-Rasterizar
r <- raster(ncol=1000, nrow=1000) #definir raster (vacío)
extent(r) <- extent(prov_Guayas) #tomar extensión espacial de capa existente
projection(r) <- proj4string(prov_Guayas) #tomar CRS de capa existente
res(r)<-c(500,500) #asignar resolucion espacial
clima_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', background=0)
#visualizar
spplot(clima_Guayas, "descripcio",main="clima vector", col.regions=rev(heat.colors(5)))
spplot(clima_ras, main="clima raster", col.regions=rev(heat.colors(20)))
spplot(hidrogeologia_Guayas, "PERMIABILI", main="Hidrogeología vector",
col.regions=rev(heat.colors(7)))
spplot(hidrogeologia_ras, main="hidrogeologia raster",
col.regions=rev(heat.colors(20)))
spplot(prov_Guayas, "DPA_DESPRO",col="white", fill="grey",
sp.layout = list("sp.polygons", zonas_inundables_Guayas,
col="blue", fill="blue"),colorkey=FALSE)
spplot(zonas_inundables_ras, main="area inundacion raster",
col.regions=c("grey", "blue"), colorkey=FALSE)
#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)spplot(algebra_mapas, col.regions=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.4650032#alternativa 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.6014852#alternativa 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?
Has llegado al final del material!