Curso de Análise de Dados Geográficos com R
ISA, Fevereiro de 2015
Módulos: SIGs com R
Interpolação espacial
Manuel Campagnolo
Instituto Superior de Agronomia, Universidade de Lisboa
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
1 / 164
Dados geográficos, ambiente R e aplicações
SIG
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
2 / 164
Processar informação geográfica em ambiente R
O objectivo do módulo é mostrar como se podem manipular dados
geográficos em R. Em particular,
1
2
Ler, associar e alterar sistemas de coordenadas de um conjunto
de dados geográficos;
Ler, criar e exportar conjuntos de dados geográficos de tipo:
1
2
matricial (“raster”), e.g. formato geotiff;
vectorial, e.g. formato shapefile.
Links sugeridos:
1
The Comprehensive R Archive Network’s task view “Analysis of
Spatial Data”, by Roger Bivand
2
Edzer Pebesma and Roger Bivand, Classes and methods for
spatial data in R, R news, 5/2 (2005) 9–14
3
Robert J. Hijmans, Introduction to the raster package, 2014
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
3 / 164
Dados geográficos
As funções base do R permitem ler, visualizar e analisar dados. Neste
módulo, o interesse é em dados geográficos, isto é, observações
associadas a localizações geográficas.
Dados geográficos têm então duas componentes:
A localização das observações, que é definida num certo CRS:
sistema de coordenadas (geográficas ou cartográficas);
Os valores associados a essas observações que podem ser
representados por variáveis contínuas (e.g. temperatura,
precipitação, elevação, índice de área foliar, área, etc) ou
categoriais (e.g. nome de freguesia, tipo de ocupação do solo,
tipo de solo, tipo de estrada, estado operacional de uma estação
meteorológica, etc)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
4 / 164
Packages R para dados geográficos
Existem muitos packages disponíveis em R para dados geográficos
(ver The Comprehensive R Archive Network’s task view “Analysis of
Spatial Data”, by Roger Bivand para uma descrição do conjunto de
packs disponíveis).
Essencialmente, os packages auxiliam o utilizador a:
1
importar para ambiente R dados geográficos e exportar dados
geográficos que resultam da análise feita em R;
2
analisar dados geográficos em R.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
5 / 164
Relação entre dados geográficos, ambiente R e aplicações SIG
Por exemplo, uma imagem Landsat é descarregada do site usgs.gov
em formato GeoTIFF. Pode ser lida em ambiente R através do
comando landsat<-raster(nomeFicheiro) que devolve um
objecto “raster”. É legítmo fazer então m<-as.matrix(landsat).
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
6 / 164
Sistema de coordenadas (CRS): descrição proj.4
O sistema de coordenadas pode ser escrito de várias formas. Em R é
comum usar-se a descrição proj.4 como nos exemplos seguintes:
1
Coordenadas geográficas WGS84 :
"+proj=longlat +ellps=WGS84 +datum=WGS84";
2
Coordenadas projectadas ETRS-PT-TM06
"+ellps=GRS80 +towgs84=0,0,0 +proj=tmerc
+lat_0=39d40’05.73’’N+lon_0=08d07’59.19’’W";
3
Coordenadas “militares” do IGeoE
"+proj=tmerc +towgs84=-304.046,-60.576,103.64,0,0,0
+k=1 +x_0=200000 +y_0=300000 +ellps=intl
+pm=lisbon +units=m".
Nota: uma boa fonte para descrições proj.4 e outras descrições de
sistemas de coordenadas é o site
http://spatialreference.org/.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
7 / 164
Estruturas para dados geográficos matriciais
(raster) em R: o package raster.
Transformações de sistemas de coordenadas
em R
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
8 / 164
Dados matriciais: leitura e escrita
A leitura de dados do tipo matricial pode ser feita com as funções:
1
raster do package raster; esta função converte dados
matriciais (no ambiente R ou fora do ambiente R, i.e. num
ficheiro) num objecto RasterLayer;
2
readGDAL do package rgdal; Esta função permite ler ficheiros
de dados matriciais e tipicamente devolve um objecto
SpatialGridDataFrame; GDALinfo devolve os parâmetros do
conjunto de dados.
A escrita de dados pode ser feita com
1
writeRaster do package raster;
2
writeGDAL do package rgdal;
Nota: raster é tipicamente mais eficiente do que readGDAL para a
leitura de conjuntos de dados de grande dimensão.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
9 / 164
Raster representation
All raster data have essentially the same structure: a regularly spaced
array of data. Each element of the array is called a grid cell and has a
numerical value or a no-value (No data, Null, ...).
Figure: Illustration of the general structure of raster data. In the simplest
case, geographic raster data are georeferenced by the coordinates of a single
corner of the array and the cell size (spatial resolution).
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
10 / 164
Raster Data Formats
Raster data consists of spatially coherent, digitally stored numerical
data, collected from sensors, scanners, or in other ways numerically
derived.
Raster data values, as read in from a file, are organized by software
into two dimensional arrays, the indices of the arrays being used as
coordinates (note: there may also be additional indices for
multispectral data).
There are many different raster formats used by different softwares
and raster image producers. Most GIS can read a multitude of raster
formats and export raster data in a variety of formats.
Currently, the most popular and versatile raster data format is Tag
Image File Format (TIFF) , with filename extension .tiff or .tif.
The TIFF file format allows for a flexible set of information fields (tags).
It consists of a single file up to 4Gb and it is a lossless format.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
11 / 164
Raster Geographical Data Models: the GeoTIFF format
The GeoTIFF data format is the TIFF format plus geolocation tags.
The geolocation is defined as follows.
1
Raster space R or image space, used to reference the pixel
values in an image: it is a regular grid with (0,0) denoting the
upper-left corner; The grid has a given number of columns and
number of rows; Each cell value can be though as associated to
an area, i.e. the first cell value is associated to the cell with center
(0.5,0.5), or to a point, i.e. the first cell value is associated to (0,0);
2
Model space M, used to reference points on the Earth: it has a
Coordinate Reference System (CRS) which may include vertical
coordinates.
3
A transformation between R and M, which can be defined as one
tie point and pixel size (in x and y ), a transformation matrix, a list
of ground control points (typically 3 corners of the image), or
rational polynomial coefficients for orthorectification.
The resolution is the size of a cell after transformation. The data type
of the data can be integer or “float” and of a given size (bits).
For details see http://trac.osgeo.org/geotiff/
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
12 / 164
Dados matriciais: leitura e escrita com package raster:
1
2
3
4
5
6
7
# c a r r e g a r package :
library ( raster )
# d e f i n i r a pasta de t r a b a l h o :
setwd ( "C: / . . . / cursoR " )
# ler ficheiro :
pan <− r a s t e r ( " landsat8pan . t i f " )
p l o t ( pan )
O objecto pan é um objecto R de classe RasterLayer:
1
2
3
4
pan # descreve o c o n j u n t o de dados
pan@crs # devolve o sistema de coordenadas em f o r m a t o p r o j . 4
que é + p r o j =utm +zone=29 +datum=WGS84 + u n i t s =m +no_ d e f s +
e l l p s =WGS84 +towgs84 =0 ,0 ,0
# e x p o r t a r para f i c h e i r o GeoTIFF
w r i t e R a s t e r ( r , f i l e = " o u t . t i f " , f o r m a t = " G T i f f " , o v e r w r i t e =TRUE)
Nota. O package raster define várias classes de objectos
(RasterLayer, Extent, ...). Cada classe contém slots que servem a
armazenar informação da classe. O acesso a slots faz-se com @ (e.g.
pan@crs) ou com a função slot.
Manuel Campagnolo
(ISA/ULisboa)
SIGs com Rraster são
5 a 10listadas
de Fevereiro de
2015
13 / 164
Nota:
as opções
de plot para objectos
com
Análise de dados com package raster: imagem pancromática
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
14 / 164
Análise de dados com package raster: valores das células
A imagem obtida com plot(pan) usa uma paleta de cores por
omissão. Podemos, como é habitual com imagens pancromáticas,
usar tons de cinzento:
1
2
3
4
5
# c r i a r v e c t o r de 100 t o n s de c i n z e n t o e n t r e p r e t o ( 0 ) e branco
(1)
cores<−gray ( seq ( 0 , 1 , l e n g t h . o u t =100) )
p l o t ( pan , c o l = cores )
# t i r a r c a i x a e c o l o c a r a legenda na h o r i z o n t a l
p l o t ( pan , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE , c o l =cores
, h o r i z o n t a l =TRUE)
A imagem está escura dado que plot transforma linearmente os
valores das células num intervalo [0,1], ou seja, transforma o valor
mínimo em 0, e o valor máximo em 1, e todos os valores intermédios
através da recta correspondente.
Os valores das células mínimo e máximo podem ser obtidos com:
1
2
pan@data@min
pan@data@max
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
15 / 164
Raster data sets: histogram and contrast enhancement
A raster dataset with grid cell values {z1 , . . . , zn } has range [zmin , zmax ].
The j-th grid cell, with value zj , is displayed with a some intensity Ij on
the screen which ranges from Imin (say, 0) to Imax (say, 1). Usually, Imin
is displayed as dark and Imax is displayed as light.
With no contrast enhancement, the intensity I is a linear function of z
over the whole [zmin , zmax ]; hence, the histogram of the I’s has the
same shape as the histogram of the z’s. Contrast enhancement is
defined by a (increasing, non-linear) function I = f (z).
The simplest case is to select a
new range [znew min , znew max ] and
define f continuous such that
f (z) = Imin , when z < znew min ,
f (z) = Imax , when z > znew max ,
and f linear in between.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
16 / 164
Análise de dados com package raster: valores das células
Os valores das células de um objecto RasterLayer podem ser lidos
para o ambiente R com a função values:
1
v<−v a l u e s ( pan ) # v e c t o r numérico
Nota: a leitura de cdg matriciais com raster não coloca a totalidade
dos dados na sessão R por uma questão de eficiência; a função
values lê os valores no ficheiro "landsat8pan.tif".
Agora é possível fazer o histograma e perceber como é que as
intensidades na imagem são calculadas:
1
2
3
4
5
6
7
o u t<− h i s t ( v , main= " v a l o r e s vs i n t e n s i d a d e " )
f q s<−o u t $ counts # f r e q u ê n c i a s a b s o l u t a s das c l a s s e s
# r e p r e s e n t a r a t r a n s f or m aç ã o l i n e a r em i n t e n s i d a d e
l i n e s ( x=range ( v ) , y=range ( f q s ) , c o l = " b l u e " )
# r e p r e s e n t a r os v a l o r e s de i n t e n s i d a d e
a x i s ( 4 , a t =seq ( min ( f q s ) ,max ( f q s ) , l e n g t h . o u t =5) , l a b e l s =seq ( 0 , 1 ,
l e n g t h . o u t =5) )
mtext ( " i n t e n s i d a d e " , s i d e =4 , l i n e = −1.5 , c o l = " b l u e " )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
17 / 164
Análise de dados com package raster: valores das células
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
18 / 164
Análise de dados com package raster: valores das células
Uma forma de aumentar o contraste é saturar a imagem (para valores
baixos, i.e. escuros) ou valores elevados. Neste caso, a observação
do histograma sugere saturar a imagem para valores elevados.
Pode usar-se o argumento zlim fixando a gama de z’s:
1
p l o t ( pan , z l i m =c ( 0 , 2 5 0 0 0 ) , x a x t = " n " , y a x t = " n " , box=FALSE , axes
=FALSE , c o l =cores , h o r i z o n t a l =TRUE)
Em alternativa pode saturar-se pelo quantil 0.999:
1
p l o t ( pan , z l i m =c ( 0 , q u a n t i l e ( v , . 9 9 9 ) ) , x a x t = " n " , y a x t = " n " , box=
FALSE , axes=FALSE , c o l =cores , legend=FALSE )
Em analogia a zlim pode usar-se xlim e ylim:
1
p l o t ( pan , x l i m =c (570000 ,575000) , y l i m =c (4315000 ,4325000) , z l i m =
c ( 0 , 2 5 0 0 0 ) , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE , c o l =
cores , legend=FALSE )
coordinates devolve as coordenadas; res devolve resolução:
1
head ( c o o r d i n a t e s ( pan ) ) ; r e s ( pan )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
19 / 164
Análise de dados com package raster: mudança de contraste
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
20 / 164
Dados matriciais: interface com “rato” e extent
O ambiente R é muito limitado para manipular interactivamente dados
geográficos. No entanto existem funções que permitem alguma
interactividade tais como zoom, click, locator ou drawExtent.
A função drawExtent permite usar o rato sobre a imagem para
definir um objecto de classe extent:
1
2
3
4
p l o t ( pan , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE , legend=
FALSE , c o l = cores )
box <− drawExtent ( ) # c l i c a r em 2 cantos
# d e f i n i r a extensão da imagem usando box
p l o t ( pan , e x t =box , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE ,
legend=FALSE , c o l = cores )
Pode manipular-se directamente o objecto da classe extent para
alterar a extensão da figura. No exemplo abaixo, usa-se só a metade
central de box:
1
p l o t ( pan , e x t =box / 2 , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE ,
legend=FALSE , c o l =cores )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
21 / 164
Dados matriciais: interface com “rato” e extent
A função zoom permite fazer um "zoom" numa região rectangular:
1
2
p l o t ( pan , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE , legend=
FALSE , c o l = cores )
zoom ( pan , new=FALSE , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE ,
legend=FALSE , c o l =cores ) # c l i c a r em 2 cantos
Pode aceder-se aos valores das células interactivamente com click:
1
c l i c k ( pan ) # devolve v a l o r de p i x e l s ; p a r a r com <ESC>
Pode digitalizar-se um polígono com locator
1
2
p t s <− l o c a t o r ( ) # devolve uma l i s t a de coords ; p a r a r com <ESC>
polygon ( c b i n d ( p t s $x , p t s $y ) , b o r d e r = " red " ) # a d i c i o n a à imagem o
p o l í g o n o formado p e l o s pontos d i g i t a l i z a d o s
Pode usar-se extent ou algum objecto que possa ser interpretado
como tal para recortar pan e criar um novo RasterLayer:
1
2
pan . crop <− crop ( x=pan , y= p t s ) # a nova extensão é dada por p t s
p l o t ( pan . crop , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE ,
legend=FALSE , c o l = cores )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
22 / 164
Análise de dados com package raster: alterar sistema de
coordenadas
A re-projecção de um cdg matricial, com um dado CRS, num novo cdg
matricial com outro CRS envolve três aspectos:
1
transformação de coordenadas f (x, y ) = (u, v );
2
características da grelha de output (extensão, resolução);
3
critério de re-amostragem.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
23 / 164
Análise de dados com package raster: alterar sistema de
coordenadas
O package raster contém a função projectRaster que permite
re-projectar os dados, definir a nova resolução e extensão, e escolher
o critério de re-amostragem.
projectRaster usa os sistemas de coordenadas (CRS) do input e
do output. Como visto atrás, a forma mais simples de definir o CRS é
através de uma descrição proj.4.
A descrição proj.4 contém em particular:
1
A definição da projecção cartográfica, e.g. +proj=tmerc;
2
A definição do datum, e.g. +datum=WGS84;
3
A descrição da transformação de datum relativamente ao datum
WGS84, e.g.
+towgs84=-304.046,-60.576,103.64,0,0,0,0 ou
+nadgrids=ptLX_e89.gsb.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
24 / 164
Análise de dados com package raster: alterar sistema de
coordenadas
Exemplo 1: re-projectar indicando apenas o CRS do output; neste
caso, a extensão e resolução são iguais às do input; o critério de
re-amostragem é bi-linear por omissão.
1
2
igeoe <− " + p r o j =tmerc + l a t _0=39.66666666666666 +towgs84
= −304.046 , −60.576 ,103.64 ,0 ,0 ,0 ,0 + l o n _0=1 +k=1 +x_0=200000
+y_0=300000 + e l l p s = i n t l +pm= l i s b o n + u n i t s =m" # Coordenadas
‘ ‘ militares ’ ’
m <− p r o j e c t R a s t e r ( pan , c r s =igeoe )
Exemplo 2: criar primeiro um objecto RasterLayer vazio com a
extensão e resolução desejados, e em seguida fazer a re-projecção do
input com projectRaster ; no exemplo é também escolhido o
critério de re-amostragem com o argumento method:
1
2
3
v a z i o <− r a s t e r ( xmn=195000 , xmx=205000 , ymn=230000 , ymx=245000 ,
r e s o l u t i o n =100 , c r s =igeoe )
m2 <− p r o j e c t R a s t e r ( from=pan , t o = vazio , method= " ngb " )
m3 <− p r o j e c t R a s t e r ( from=pan , t o = vazio , method= " b i l i n e a r " )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
25 / 164
Análise de dados com package raster: alterar sistema de
coordenadas
Ilustração: comparar input e output da re-projecção.
Representar dados originais em coordenadas UTM, zona 29:
1
2
3
4
5
par ( mfrow=c ( 2 , 1 ) )
p l o t ( pan , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE , legend=
FALSE , c o l =cores , z l i m =c ( 0 , 2 5 0 0 0 ) )
xy<−as . v e c t o r ( pan@extent )
# xpd para poder e s c r e v e r f o r a da imagem
t e x t ( xy [ c ( 1 , 2 , 1 , 1 ) ] , xy [ c ( 3 , 3 , 3 , 4 ) ] , round ( xy ) , pos=c ( 1 , 1 , 2 , 2 ) , xpd
=TRUE)
Representar a imagem transformada para coordenadas militares, com
uma transformação de datum definida por três parâmetros (igeoe):
1
2
3
p l o t (m, x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE , legend=FALSE
, c o l =cores , z l i m =c ( 0 , 2 5 0 0 0 ) )
xy<−as . v e c t o r ( m@extent ) # extensão xmin , xmax , ymin , ymax
t e x t ( xy [ c ( 1 , 2 , 1 , 1 ) ] , xy [ c ( 3 , 3 , 3 , 4 ) ] , round ( xy ) , pos=c ( 1 , 1 , 2 , 2 ) , xpd
=TRUE)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
26 / 164
Análise de dados com package raster: sistemas de coordenadas
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
27 / 164
Análise de dados com package raster: alterar sistema de
coordenadas
A transformação de datum (que é sempre relativa a WGS84) pode ser
feita de duas formas:
1
através de um conjunto de parâmetros (transformação 3D
conhecida por Helmert ou Bursa-Wolf), e.g.
+towgs84=-304.046,-60.576,103.64,0,0,0,0;
2
através de grelhas (conhecidas por "NAD grids") , e.g.
+nadgrids=ptLX_e89.gsb.
Neste caso, é necessário:
1
2
3
colocar o ficheiro com a informação necessária (e.g.
"ptLX_e89.gsb") na pasta indicada por
system.file("proj", package = "rgdal");
indicar na definição do CRS que a transformação de datum é do
tipo "nadgrid" com o parâmetro +nadgrids=ptLX_e89.gsb;
aplicar a re-projecção da forma habitual.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
28 / 164
Análise de dados com package raster: projecção com grelhas
Exemplo: Re-projectar cdg pan pelo método das grelhas.
Definir CRS com parâmetro +nadgrids, e re-projectar:
1
2
igeoe . g r i d <− CRS( " + p r o j =tmerc + l a t _0=39.66666666666666 +
n a d g r i d s =ptLX_e89 . gsb + l o n _0=1 +k=1 +x_0=200000 +y_0=300000
+ e l l p s = i n t l +pm= l i s b o n + u n i t s =m" )
m. g r i d <− p r o j e c t R a s t e r ( pan , c r s =igeoe . g r i d )
Representar o resultado:
1
2
3
p l o t (m. g r i d , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE , legend=
FALSE , c o l =cores , z l i m =c ( 0 , 2 5 0 0 0 ) )
xy<−as . v e c t o r (m. gr id @extent )
t e x t ( xy [ c ( 1 , 2 , 1 , 1 ) ] , xy [ c ( 3 , 3 , 3 , 4 ) ] , round ( xy ) , pos=c ( 1 , 1 , 2 , 2 ) , xpd
=TRUE)
Sugestão: pode comparar com o cdg m, obtido através do CRS
igeoe, que usa 3 parâmetros para a transformação de datum.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
29 / 164
Análise de dados com package raster: projecção com grelhas
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
30 / 164
Análise de dados com package raster: Exemplo de utilização de
projectRaster para determinar a localização de um ponto num
outro CRS
Usando o objecto pan, em coordenadas UTM zona 29, que
representa uma região em redor da albufeira da Barragem de
Montargil, determine:
1
a gama de latitudes e longitudes que corresponde a essa
imagem;
2
a localização da parede da barragem em lat/long.
Sugestão: Use projectRaster e o CRS de coordenadas
geográficas lat/long para obter as novas coordenadas. Para localizar
bem a parede da barragem use zoom e para identificar use locator.
Solução: aproximadamente lat=39.05347 N e long=8.176397 W.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
31 / 164
Análise de dados com package raster: Exercício
Usando o objecto pan, que representa uma região em redor da
albufeira da Barragem de Montargil, determine:
1
a gama de latitudes e longitudes que corresponde a essa
imagem;
2
a localização da parede da barragem em lat/long.
Sugestão: Use projectRaster e o CRS de coordenadas
geográficas lat/long para obter as novas coordenadas. Para localizar
bem a parede da barragem use zoom e para identificar use locator.
Solução: aproximadamente lat=39.05347 N e long=8.176397 W.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
32 / 164
Análise de dados com package raster: Exercício
Neste exercício pretende-se comparar os valores de elevação com os
valores de índice de vegetação “NDVI” para a área protegida
Sintra-Cascais. Os dados de input são:
1
Dados NDVI derivados de uma imagem Landsat 8 de Março de
2014 com CRS UTM zona 29 e resolução de 30m – ficheiro
“ndviSintra.tif”;
2
Modelo digital de elevações (MDE) SRTM com CRS WGS84 e
resolução de 3-arc segundos – ficheiro “n38_w010_3arc_v2.tif”;
3
Limite da área protegida em CRS ETRS-PT-TM06, extraído de
cartas disponíveis na página do ICNF – ficheiro
“limite.AP.SintraCascais.ETRS.txt”.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
33 / 164
Análise de dados matriciais: exercício AP Sintra-Cascais
Ler ficheiro com as coordenadas (ETRS) do limite da área protegida
para matriz:
1
l i m i t e . ap <− as . m a t r i x ( read . t a b l e ( f i l e = " l i m i t e . AP . S i n t r a C a s c a i s
. ETRS . t x t " , header=FALSE ) )
Ler ficheiro com dados de elevação (MDE):
1
mde <− r a s t e r ( " n38_w010_3 a r c _v2 . t i f " )
Converter MDE para ETRS e recortar pelo limite da AP:
1
2
3
e t r s <− " + p r o j =tmerc + l a t _0=39.6682583 + l o n _0=−8.1331083 +k=1 +
x_0=0 +y_0=0 + e l l p s =GRS80 + u n i t s =m"
mde . e t r s <− p r o j e c t R a s t e r ( mde , c r s = e t r s )
mde . e t r s <− crop (mde . e t r s , l i m i t e . ap )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
34 / 164
Análise de dados matriciais: exercício AP Sintra-Cascais
Ler imagem de índice de vegetação NDVI:
1
n d v i<− r a s t e r ( " n d v i S i n t r a . t i f " )
Re-projectar NDVI para ETRS alinhando a imagem com mde.etrs. O
alinhamento é feito escolhendo como extensão e como resolução do
novo RasterLayer a extensão e a resolução de mde.etrs:
1
2
3
n d v i<− r a s t e r ( " n d v i S i n t r a . t i f " )
v a z i o<− r a s t e r ( e x t =mde . etrs@extent , c r s =mde . etrs@crs , r e s o l u t i o n =
r e s (mde . e t r s ) )
n d v i . e t r s<−p r o j e c t R a s t e r ( from= ndvi , t o = v a z i o )
Nota: para re-projectar um cdg sobre a quadrícula de outro cdg
matricial, não é necessário conhecer a extensão, resolução e CRS
desse cdg: basta extrair essa informação com @extent, @crs, e com
a função res.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
35 / 164
Análise de dados matriciais: exercício AP Sintra-Cascais
Construir 2 transectos perpendiculares que passam pelo centro da AP
e pelos eixos principais correspondentes a essa área.
Determinar os eixos principais, dados pelos vectores próprios da
matrix de variância-covariância das coordenadas do limite da
área:
1
1
2
cp1<−eigen ( cov ( l i m i t e . ap ) ) $ v e c t o r s [ , 1 ] # 1a componente
principal
cp2<−eigen ( cov ( l i m i t e . ap ) ) $ v e c t o r s [ , 2 ] # 2a componente
principal
Extrair uma amostra de pontos ao longo desses eixos:
2
1
2
3
4
xy0<−a p p l y ( l i m i t e . ap , 2 , mean ) # ponto c e n t r a l da AP
xs<−unique ( c o o r d i n a t e s (mde . e t r s ) [ , 1 ] ) # coord x das c é l u l a s
amostra<−c b i n d ( xs , xy0 [ 2 ] + cp1 [ 2 ] / cp1 [ 1 ] * ( xs−xy0 [ 1 ] ) ) #
coord x e y
amostra<−r b i n d ( amostra , c b i n d ( xs , xy0 [ 2 ] + cp2 [ 2 ] / cp2 [ 1 ] * ( xs−
xy0 [ 1 ] ) ) )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
36 / 164
Análise de dados matriciais: exercício AP Sintra-Cascais
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
37 / 164
Análise de dados matriciais: exercício AP Sintra-Cascais
Extrair os valores dos cdg nmatriciais com a função extract,
aplicada à amostra criada anteriormente.
1
2
y<− e x t r a c t ( n d v i . e t r s , amostra )
x<− e x t r a c t (mde . e t r s , amostra )
Pode agora construir-se um gráfico e ajustar uma curva y = a x b aos
dados:
1
2
3
4
5
p l o t ( y~x , x l a b = " elevação (m) " , y l a b = " n d v i " , main= "AP S i n t r a −
Cascais " , cex . l a b = 1 . 3 )
xx<−l o g ( x [ x > 0 ] ) # t r a n f or m a çã o l o g a r í t m i c a
yy<−l o g ( y [ x > 0 ] )
a j u s t<−lm ( yy~xx ) # a j u s t a r regressão
curve ( exp ( a j u s t $ c o e f [ 1 ] ) * x ^ a j u s t $ c o e f [ 2 ] , add=TRUE, c o l = " red " )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
38 / 164
Análise de dados matriciais: exercício AP Sintra-Cascais
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
39 / 164
Processamento de composições de imagems
com o package raster
Filtros de janela móvel
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
40 / 164
Processamento de composições de imagems com o package raster
e filtros
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
41 / 164
Análise de dados com package raster: composição de imagens
A classe raster contém objectos RasterBrick e RasterStack
que suportam múltiplas imagens. Cada imagem deve ter a mesma
extensão e resolução. RasterStack é mais flexível (por exemplo,
pode ser formado lendo vários ficheiros) e RasterBrick é uma
estrutura de dados mais rígida mas que permite um processamento
mais eficiente.
No próximo exemplo usa-se a função stack para ler várias bandas
Landsat da mesma região (cada banda é um ficheiro GeoTIFF).
Depois, converte-se o resultado para um RasterBrick . A função
plotRGB permite escolher as bandas da composição colorida.
1
2
3
4
5
6
7
8
f i c h s <− l i s t . f i l e s ( p a t t e r n = " banda " )
s <− s t a c k ( as . l i s t ( f i c h s ) )
s # devolve sumário dos dados
b <− b r i c k ( s )
b # idem
n l a y e r s ( b ) # devolve número de camadas em b
box <−e x t e n t ( c (500000 , 520000 , 4310000 , 4330000) )
plotRGB ( b , r =4 ,g=3 ,b=2 , s t r e t c h = " l i n " , e x t =box ) #RGB=432
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
42 / 164
Análise de dados com package raster: Composição Landsat
RGB=432, Ribatejo
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
43 / 164
Análise de dados com package raster: dados Google maps
É possível, em ambiente R ter acesso a várias fontes de dados
matriciais, e em particular, a imagens Google Maps. Existem vários
packages que permitem o acesso e "download" dessas imagens, que
ficam disponíveis na sessão R e podem ser manipuladas como
qualquer outro conjunto de dados do R. Em particular, pode aceder-se
a cada um dos canais RGB dessas imagens.
Um dos packages com essa funcionalidade é o package dismo, que
pode ser instalado e carregado com:
1
2
i n s t a l l . packages ( " dismo " , repos= " h t t p : / / R−Forge . R−p r o j e c t . org " )
l i b r a r y ( dismo )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
44 / 164
Análise de dados com package raster: dados Google maps
Vamos escolher uma pequena zona (por forma a ter melhor detalhe)
para a qual se pretende ter as imagens de muito alta resolução Google
Maps, na parte central do stack b (no Ribatejo) usado anteriormente.
Em primeiro vamos verificar que b tem coordenadas UTM, zona 29:
1
b@crs
Vamos extrair uma pequena região da imagem Landsat – e poderia ter
sido obtido com drawExtent():
1
2
e<−e x t e n t ( as . v e c t o r ( c ( 513316 , 514509 , 4317568 ,4318417) ) )
l a n d s a t<−crop ( b , e )
A imagem Google é extraída com coordenadas lat/long:
1
2
3
4
5
wgs84<−" + p r o j = l o n g l a t + e l l p s =WGS84 +datum=WGS84"
aux<−p r o j e c t R a s t e r ( l a n d s a t , c r s =wgs84 )
# o s l o t @extent devolve um o b j e c t o de c l a s s e e x t e n t
aux@extent # coordenadas l a t / l o n g do o b j e c t o l a n d s a t
zona . c e n t r a l<−aux@extent / 2 # coordenadas l a t / l o n g da p a r t e
c e n t r a l da extensão dada por aux@extent
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
45 / 164
Análise de dados com package raster: dados Google maps
Sabendo a extensão da zona em que queremos extrair a imagem
Google, pode usar-se então a função gmap do package dismo para
aceder a essa imagem. O argumento rgb permite extrair a
composição colorida (rgb=FALSE devolve um RasterLayer) ou as
três bandas (rgb=TRUE devolve um RasterBrick). O argumento
scale pode ser 1 ou 2 (resolução duas vezes melhor).
1
gm. wgs<−gmap ( zona . c e n t r a l , t y p e = " s a t e l l i t e " , l o n l a t =TRUE, rgb=TRUE
, s c a l e =2)
Executar
1
gm. wgs
para verificar que o CRS de gm.wgs é WGS84, e que as bandas se
designam por "red", "green", "blue".
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
46 / 164
Análise de dados com package raster: dados Google maps
Em seguida vamos verificar visualmente que as imagens gm.wgs
(CRS WGS84) e landsat (CRS UTM, zona 29) se sobrepõem
correctamente.
1
Re-projectar os dados num novo CRS: ETRS-PT-TM06.
1
2
3
e t r s<−" + p r o j =tmerc + l a t _0=39.6682583 + l o n _0=−8.1331083 +k=1
+x_0=0 +y_0=0 + e l l p s =GRS80 + u n i t s =m"
gm. e t r s<−p r o j e c t R a s t e r (gm. wgs , c r s = e t r s )
l a n d s a t . e t r s<−p r o j e c t R a s t e r ( l a n d s a t , c r s = e t r s )
Em seguida, construir a sobreposição usando parâmetro alpha
que varia entre 0, para transparente, e 255, para opaco. A
composição landsat é RGB=432
2
1
2
3
4
par ( mfrow=c ( 1 , 1 ) )
plotRGB (gm. e t r s , r =3 ,g=2 ,b=1 , s t r e t c h = " l i n " ) #RGB
# composição c o l o r i d a Landsat RGB=432
plotRGB ( l a n d s a t . e t r s , r =4 ,g=3 ,b=2 , s t r e t c h = " l i n " , e x t =gm.
etrs@extent , alpha =100 ,add=TRUE)
Sugestão: repetir o procedimento para RGB=742.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
47 / 164
Análise de dados com package raster: Sobreposição de uma
composição Landsat RGB=432 sobre um mapa de alta resolução
(Google maps) no Ribatejo
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
48 / 164
Análise de dados com package raster: álgebra de imagens
Os valores de um objecto raster podem ser acedidos com a função
values que devolve um vector numérico. Mas é normalmente mais
simples usar directamente o objecto raster. No exemplo abaixo
calcula-se o NDVI a partir das bandas 3 e 4 de RasterBrick b.
1
n d v i <− ( b$banda4 − b$banda3 ) / ( b$banda4 + b$banda3 )
Os valores são manipulados como numa matriz, e.g. determinar a
proporção de pixels de ndvi que têm valor NA, ou alterar valores:
1
2
n c e l l ( n d v i [ i s . na ( n d v i ) ] ) / n c e l l ( n d v i ) #pode usar l e n g t h em vez
de n c e l l
n d v i [ ndvi <0] <− 0 # s u b s t i t u i ç ã o de v a l o r
Existem funções de mais alto nível tais como reclassify, subs ou
cut que permitem alterar os valores das células usando uma "lookup
table". Abaixo reclassifica-se ndvi de acordo com a transformação
] − 1, 0] → 1, ]0, .5] → 2 e ].5, 1] → 3:
1
2
t a b e l a <− c b i n d ( c ( − 1 , . 2 , . 5 ) , c ( . 2 , . 5 , 1 ) , 1 : 3 )
n d v i . c <− r e c l a s s i f y ( ndvi , r c l = t a b e l a , r i g h t =TRUE)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
49 / 164
Análise de dados com package raster: álgebra de imagens
Um conjunto de funções aplica-se a RasterBrick e RasterStack
célula a célula devolvendo um novo raster
1
2
min ( b ) # devolve o mínimo nas 6 bandas
range ( b ) # devolve R a s t e r B r i c k com min e max
Para obter uma estatística para cada banda usa-se cellStats:
1
c e l l S t a t s ( b , " mean " ) # devolve um v e c t o r
Operações numéricas ou lógicas sobre um ou mais objectos raster
podem ser realizadas com as funções calc e overlay
respectivamente; fun tem que aceitar o mesmo número de
argumentos que o número de camadas de RasterBrick ou
RasterStack:
1
2
median ( b ) #dá e r r o
b . mediana<−o v e r l a y ( b , f u n =median , rm . na=TRUE) # l e n t o , mas
funciona
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
50 / 164
Análise de dados com package raster: mosaico de imagens
Podemos criar um mosaico a partir de um conjunto de imagens
alinhadas, com a mesma resolução e CRS. Existem para tal duas
funções no package raster:
1
a função merge que usa a primeira imagem do conjunto quando
há sobreposição;
2
a função mosaic que é mais flexível, e que permite usar qualquer
função das várias imagens quando há sobreposição.
A função tem um argumento fun que é aplicado ao conjunto de
valores de células com a mesma localização. A função indicada
por fun tem que ter obrigatoriamente um argumento na.rm para
poder ser aplicada a células com valor NA. Por exemplo,
fun=min,
fun=mean,
fun=function(x,...) quantile(x,.95).
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
51 / 164
Análise de dados com package raster: mosaico de imagens
Vamos usar dados de altimetria SRTM disponíveis on-line. Primeiro
vamos criar uma string com a URL onde estão os dados:
1
HTTP<−" h t t p : / / dds . c r . usgs . gov / srtm / v e r s i o n 2 _1 /SRTM3/ E u r a s i a /
N37W008 . h g t . z i p " # cobre SE de P o r t u g a l
Depois, usa-se a função do R download.file para descarregar o
ficheiro para a pasta de trabalho:
1
download . f i l e ( u r l =HTTP, d e s t f i l e = "N37W008 . h g t . z i p " ,mode= "wb" )
Como o ficheiro está comprimido, usa-se unzip para descomprimir
na pasta de trabalho:
1
u n z i p ( z i p f i l e = "N37W008 . h g t . z i p " )
Finalmente, usa-se raster, que lê em particular o formato "hgt":
1
srtm8<− r a s t e r ( "N37W008 . h g t " )
Fazer o mesmo para descarregar o ficheiro N37W009, e obter srtm9.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
52 / 164
Análise de dados com package raster: mosaico de imagens
Para os dados do exemplo que está a ser considerado, não há
sobreposição das imagens e por isso o mosaico pode ser
simplemente obtido com a função merge:
1
srtm <− merge ( srtm8 , srtm9 )
O resultado pode ser representado com a figura en tons de cinzento
obtida com
1
p l o t ( srtm , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE , legend=
FALSE , c o l = cores )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
53 / 164
Análise de dados com package raster: Modelo digital de
elevações, Algarve
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
54 / 164
Análise de dados com package raster: mosaico de imagens
Para tornar a figura mais apelativa, podemos usar uma escala de
cores apropriada para a altimetria. O package RColorBrewer
permite definir algumas escalas de cores interessantes:
1
2
l i b r a r y ( RColorBrewer ) # para c r i a r p a l e t a s
d i s p l a y . brewer . a l l ( ) # v e r p o s s í v e i s p a l e t a s de cores
Vamos usar uma dessas paletas de cores para adicionar a imagem de
altimetria sobre a área continental. Primeiro, vamos atribuir valor NA
às células com elevação não positiva, e depois adicionar sobre a
figura anterior a carta hipsométrica com 12 classes:
1
2
srtm [ srtm <=0]<−NA
p l o t ( srtm , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE ,
h o r i z o n t a l =TRUE, c o l =brewer . p a l ( n=12 , " S p e c t r a l " ) , add=TRUE)
Adicionar a localização e nome da cidade de Lagos à imagem:
1
2
p o i n t s ( y =37.095753 , x = −8.683319 , c o l = " w h i t e " )
t e x t ( y =37.095753 , x = −8.683319 , " Lagos " , c o l = " w h i t e " , pos =3)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
55 / 164
Análise de dados com package raster: Carta hipsométrica,
Algarve
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
56 / 164
Raster datasets: digital elevation models and terrain analysis
Most GIS include operations to derive variables like slope, aspect and
flow direction from a digital elevation model in raster format.
The main idea is to estimate for each cell center (x, y ) the plane
tangent to the surface with slope Sx along the x direction and slope Sy
along the y direction. Using 8 neighbors, the usual estimates are:
Sx =
(z3 −z1 ) + 2 (z6 −z4 ) + (z9 −z7 )
(z1 −z7 ) + 2 (z2 −z8 ) + (z3 −z9 )
, Sy =
.
8 × resolution along x
8 × resolution along y
For this example where the resolution
is 20m the estimates are:
Manuel Campagnolo (ISA/ULisboa)
1
Sx = (1 + 2 × 3 + 3)/(8 × 20) =
0.0625, and
2
Sy = (39+2×35+37)/(8×20) =
0.9125.
SIGs com R
5 a 10 de Fevereiro de 2015
57 / 164
Raster datasets: digital elevation models and terrain analysis
(cont.)
The following variables are derived from the estimates Sx and Sy .
1
The flow direction, i.e. the direction of steepest descent, given by
the vector (−Sx , −Sy );
q
2
The slope is the slope along the flow direction S = Sx2 + Sy2 ;
3
The aspect is the flow direction expressed as the azimuth (relative
to the North, in degrees) or a cardinal (N,E,S,W) or intercardinal
(NE, SE, SW, NW) direction. The azimuth is arctan(Sx /Sy ) plus
0o (1st Quadrant), 180o (2nd or 3rd Quadrants), or 360o (4th
Quadrant). If Sx = Sy = 0 the aspect is undefined (flat terrain).
azimuth
337.5o – 22.5o
67.5o – 112.5o
157.5o – 202.5o
247.5o – 292.5o
Manuel Campagnolo (ISA/ULisboa)
direction
N
E
S
W
azimuth
22.5o – 67.5o
112.5o – 157.5o
202.5o – 247.5o
292.5o – 337.5o
SIGs com R
direction
NE
SE
SW
NW
5 a 10 de Fevereiro de 2015
58 / 164
Raster datasets: digital elevation models and terrain analysis
(cont.)
In this example,
1
2
3
In this example,
Sx =0.0625, Sy =0.9125,
√
S= 0.06252 +0.91252 =
0.9146379,
1
2
Since arctan(Sx /Sy )=3.918o ,
the aspect is 183.918o , i.e. S.
Manuel Campagnolo (ISA/ULisboa)
3
SIGs com R
Sx =−0.36875, Sy =−0.73125,
√
S= 0.368752 +0.731252 =
0.8189647,
Since arctan(Sx /Sy )=26.76o ,
the aspect is 26.76o , i.e. NE.
5 a 10 de Fevereiro de 2015
59 / 164
Raster datasets: digital elevation models and terrain analysis
(cont.)
In addition to the previous variables, it is possible to derive hillshading
and visibility from a DEM.
Hillshading is proportional to the estimated cosine of the angle of
incidence of light over the surface (θi ). This depends of the position of
the light source (typically the Sun) which is determined by the Solar
Zenith angle (θsun ) and the Solar Azimuth angle (φsun ). Given the
aspect (φaspect ) and the slope (θslope ) in degrees,
cos θi = cos θsun cos θslope + sin θsun sin θslope cos(φsun − φaspect ).
The maximum value of hillshading corresponds to θi = 0, i.e. when the
Sun radiation is normal to surface at the location.
Visibility is a binary variable that indicates if a position (grid cell) is or
not visible from a location with coordinates (x, y , Z ). Usually the
location is indicated as a position (x, y ) and an elevation above ground
h, with Z = h + z(x, y ), where z(x, y ) is the elevation given by the
DEM for the grid cell that contains (x, y ).
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
60 / 164
Análise de dados com package raster: altimetria
O package raster contém a função terrain que permite derivar de
um modelo digital de elevações um conjunto de variáveis tais como
relevo, orientação, direcção de escoamento, índices de rugosidade do
terreno, etc.
Também contém a função hillShade para derivar a iluminação sobre
o terreno a partir de um MDE e da posição da fonte de iluminação.
Por exemplo, a inclinação em radianos pode ser obtida com o seguinte
comando (pode ser aplicado a um MDT em coordenadas geográficas
desde que a elevação seja em metros):
1
i n c<− t e r r a i n ( srtm , o p t = " s l o p e " , u n i t s = " r a d i a n s " , n e i g h b o r s =8)
Os declives (%, tangente da inclinação) são dados por:
1
2
d e c l i v e s<−100 * t a n ( v a l u e s ( i n c ) ) # em %
range ( d e c l i v e s , na . rm=TRUE) # devolve minimo e máximo
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
61 / 164
Análise de dados com package raster: Exercício
Determine a localização do ponto com declive máximo. Indique a
localização desse ponto na imagem de srtm.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
62 / 164
Análise de dados com package raster: Exercício
Determine a localização do ponto com declive máximo. Indique a
localização desse ponto na imagem de srtm.
Determinar localização do declive máximo usando as coordenadas
das células:
1
2
3
aux<−c b i n d ( c o o r d i n a t e s ( srtm ) , d e c l i v e s )
colnames ( aux )<−c ( " l o n g " , " l a t " , " d e c l i v e " )
l o c . max<−aux [ which ( aux [ , 3 ] = = max ( d e c l i v e s , na . rm=TRUE) ) , ]
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
62 / 164
Análise de dados com package raster: Exercício
Determine a localização do ponto com declive máximo. Indique a
localização desse ponto na imagem de srtm.
Determinar localização do declive máximo usando as coordenadas
das células:
1
2
3
aux<−c b i n d ( c o o r d i n a t e s ( srtm ) , d e c l i v e s )
colnames ( aux )<−c ( " l o n g " , " l a t " , " d e c l i v e " )
l o c . max<−aux [ which ( aux [ , 3 ] = = max ( d e c l i v e s , na . rm=TRUE) ) , ]
Construir imagem e localizar ponto:
1
2
p l o t ( srtm , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE , legend=
FALSE , c o l = cores )
p o i n t s ( x= l o c . max [ " l o n g " ] , y= l o c . max [ " l a t " ] , c o l = " red " , cex =2)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
62 / 164
Análise de dados com package raster: localização do ponto de
maior declive
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
63 / 164
Raster datasets: image filters
A linear filter is a standard focal image processing technique: a linear
filter replaces each grid cell value by a weighted sum of its neighbors.
The kernel defines the neighborhood and the weights.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
64 / 164
Raster datasets: linear filters
Linear filters can be designed to address problems in image
processing like the following.
1
Image smoothing;
2
Computation of image derivatives and gradiants;
3
Approximate (discretize) continuous filters (e.g. Gaussian filters);
4
Template matching, i.e. determine parts of the image that are
similar to a given template.
In Slide 57 the images Sx and Sy were derived from a DEM.
Sx and Sy can be defined with the following 3 × 3 linear filters:
1
− 8×20
0

2
SX : 
 − 8×20
1
− 8×20
0

0
Manuel Campagnolo (ISA/ULisboa)
1
8×20
2
8×20
1
8×20






SY : 

SIGs com R
1
8×20
2
8×20
1
8×20

0
0
0
1
− 8×20
2
− 8×20
1
− 8×20



5 a 10 de Fevereiro de 2015
65 / 164
Raster datasets: image filters (cont.)
Filters can be used for many other operations on rasters. When an
image exhibits “salt-and-pepper noise”, different kinds of filters can be
used to reduce that noise. Below, the image on the right is “filtered”.
In particular the following non-linear filters are used to that end:
1
the median filter, which replaces the central grid cell value by the
median of its neighbors;
2
the majority filter, which replaces the central grid cell value by the
most frequent value among its neighbors.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
66 / 164
Análise de dados com package raster: filtros não lineares
Apesar da função terrain poder ser conveniente, todas as opções
(slope, aspect, TPI, TRI, rough, etc) podem ser obtidos pela
aplicação de filtros à imagem, como descrito abaixo:
1
2
3
4
5
6
7
f <− m a t r i x ( 1 , nrow =3 , n c o l =3)
# í n d i c e de rugosidade do t e r r e n o
TRI <− f o c a l ( x , w= f , f u n = f u n c t i o n ( x , . . . ) sum ( abs ( x[−5]−x [ 5 ] ) ) /
8 , pad=TRUE, padValue=NA)
# í n d i c e de posição t o p o g r á f i c a
TPI <− f o c a l ( x , w= f , f u n = f u n c t i o n ( x , . . . ) x [ 5 ] − mean ( x [ − 5 ] ) ,
pad=TRUE, padValue=NA)
# a m p l i t u d e das d i f e r e n ç a s de elevação
rough <− f o c a l ( x , w= f , f u n = f u n c t i o n ( x , . . . ) max ( x ) − min ( x ) ,
pad=TRUE, padValue=NA, na . rm=TRUE)
Obviamente, o ambiente R permite alterar a dimensão e a expressão
do filtro.
Exercício: construir um filtro de moda para substituir o valor da célula
pelo valor mais frequente entre os vizinhos. Solução: ver script R.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
67 / 164
Análise de dados com package raster: filtros lineares
Para derivar o declive de um MDE usando filtros lineares:
1
2
3
1
2
3
4
5
Em primeiro é necessário re-projectar os dados para as
distâncias horizontais serem em metros.
A diferença de cotas deve ser dividida pela resolução na direcção
respectiva para obter a estimativa do declive;
O declive em cada direcção pode ser definido com os filtros
abaixo, em que Sx é o declive estimado Sx na direcção W-E, e Sy
é o declive estimado Sy na direcção S-N como no Slide 57.
srtm . crop<−crop ( srtm , e x t e n t ( c ( − 8 . 7 , − 8 . 4 , 3 7 . 1 , 3 7 . 4 ) ) ) #
Monchique
srtm . e t r s<− p r o j e c t R a s t e r ( srtm . crop , c r s = e t r s ) # para o b t e r
coordenadas em metros
Sx <− f o c a l ( srtm . e t r s , w= m a t r i x ( c ( −1 , −2 , −1 ,0 ,0 ,0 ,1 ,2 ,1) / ( 8 * r e s (
srtm . e t r s ) [ 1 ] ) , n c o l =3) )
Sy <− f o c a l ( srtm . e t r s , w= m a t r i x ( c (1 ,0 , −1 ,2 ,0 , −2 ,1 ,0 , −1) / ( 8 * r e s (
srtm . e t r s ) [ 2 ] ) , n c o l =3) )
d e c l i v e <− s q r t ( Sx^2+Sy ^ 2 )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
68 / 164
Análise de dados com package raster: declives na Serra de
Monchique
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
69 / 164
Conjuntos de dados geográficos vectoriais em
R: estruturas de dados do package sp:
classes de pontos, linhas e polígonos
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
70 / 164
Dados vectorias: introdução
Em SIGs representam-se tipicamente objectos no espaço geográfico
de forma vectorial (feature class). Esquematicamente,
“Feature class” = objectos geométricos + tabela de atributos
Uma “feature” corresponde a uma linha da tabela.
Source: http://gisatbrown.typepad.com/gis/files/spatialdatafiles.pdf
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
71 / 164
Dados vectoriais em R: package sp
O package sp disponibiliza classes e métodos para dados espaciais
em R. Um dos objectivos principais da criação deste package (em
2005) foi fornecer uma interface uniforme para dados espaciais que
pudesse ser usada por outros packages de análise espacial.
O package sp disponibiliza classes para dados vectoriais e matriciais.
Em seguida vão ser discutidas as classes de objectos vectoriais:
tipo
points
line
lines
rings
classe
SpatialPoints
SpatialPointsDataFrame
Line
Lines
SpatialLines
SpatialLinesDataFrame
Polygon
Polygons
SpatialPolygons
SpatialPolygonsDataFrame
Manuel Campagnolo (ISA/ULisboa)
atributos
data.frame
data.frame
data.frame
SIGs com R
contém
Spatial*
SpatialPoints*
Line list
Spatial*, Lines list
SpatialLines*
Line *
Polygon list
Spatial*, Polygons list
SpatialPolygons
5 a 10 de Fevereiro de 2015
72 / 164
Dados vectoriais em R: criar SpatialPointsDataFrame
Carregar package sp
1
2
l i b r a r y ( rgdal )
l i b r a r y ( sp )
Vamos usar os dados anteriores sobre a área protegida
Sintra-Cascais, com o CRS ETRS-PT-TM06:
1
2
3
4
5
6
l i m i t e . ap<−as . m a t r i x ( read . t a b l e ( f i l e = " l i m i t e . AP . S i n t r a C a s c a i s .
ETRS . t x t " , header=FALSE ) )
f i c h . mde<−" n38_w010_3 a r c _v2 . t i f "
mde<− r a s t e r ( f i c h . mde)
e t r s<−" + p r o j =tmerc + l a t _0=39.6682583 + l o n _0=−8.1331083 +k=1 +x_
0=0 +y_0=0 + e l l p s =GRS80 + u n i t s =m"
mde . e t r s<−p r o j e c t R a s t e r ( mde , c r s = e t r s )
mde . e t r s<−crop (mde . e t r s , l i m i t e . ap )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
73 / 164
Dados vectoriais em R: criar SpatialPointsDataFrame
O objecivo deste primeiro exercício é criar um objecto de classe
SpatialPointsDataFrame que represente a localização, a
designação, e o código postal de algumas povoações, usando o CRS
ETRS-PT-TM06.
A informação de partida é a localização em lat/long, como definido na
tabela (data.frame) seguinte:
1
2
3
4
5
6
pov<−data . frame (
designa=c ( " S i n t r a " , " Colares " , " Alcabideche " , " Almoçageme " , "
Azoia " ) ,
l o n g =c ( −9.3816589 , −9.4426866999 , −9.4098824999 ,
−9.470629899 , −9.476644970130) ,
l a t =c (38.8028687 , 38.8066307 , 38.7331569 , 38.7967822 ,
38.7697973006) ,
c o d P os t a l =c ( " 2710 " , " 2705 " , " 2646 " , " 2705 " , " 2705 " )
)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
74 / 164
Dados vectoriais em R: criar SpatialPointsDataFrame
Dado que as localizações disponíveis são em lat/long começamos por
criar um objecto SpatialPointsDataFrame no CRS WGS84:
1
2
wgs84<−" + p r o j = l o n g l a t + e l l p s =WGS84 +datum=WGS84"
xy<−c b i n d ( pov$ long , pov$ l a t ) # m a t r i x de coordenadas l o n g / l a t
Criar objecto de classe SpatialPoints; a função CRS cria um
objecto de classe CRS:
1
xy . sp <− S p a t i a l P o i n t s ( xy , p r o j 4 s t r i n g =CRS( wgs84 ) )
Criar objecto de classe SpatialPointsDataFrame:
1
pov . wgs<− S p a t i a l P o i n t s D a t a F r a m e ( xy . sp , pov [ , c ( " designa " , "
codPostal " ) ] )
Esse objecto pode agora ser re-projectado para o CRS
ETRS-PT-TM06 com a função spTransform do package sp:
1
2
e t r s<−" + p r o j =tmerc + l a t _0=39.6682583 + l o n _0=−8.1331083 +k=1 +x_
0=0 +y_0=0 + e l l p s =GRS80 + u n i t s =m"
pov . e t r s<−spTransform ( pov . wgs ,CRS( e t r s ) )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
75 / 164
Dados vectoriais em R: manipular SpatialPointsDataFrame
Pode aceder-se à tabela de atributos de um objecto
SpatialPointsDataFrame com o slot @data:
1
pov . etrs@data # devolve data . frame
Para fazer uma selecção por atributos faz-se simplesmente uma
selecção de linhas da tabela de atributos como nos exemplos abaixo:
1
pov . e t r s [ 2 , ] # devolve S p a t i a l P o i n t s D a t a F r a m e com uma ú n i c a
povoação ( a 2a : Colares )
Naturalmente pode fazer-se uma selecção por atributos com uma
expressão lógica; como as colunas da tabela são do tipo "factor",
primeiro converte-se codPostal em numerico:
1
cps<−as . numeric ( as . c h a r a c t e r ( pov . etrs@data$ c o d P o s ta l ) )
e finalmente faz-se a selecção como se fosse uma tabela:
1
pov . e t r s [ cps >2700 ,] # devolve S p a t i a l P o i n t s D a t a F r a m e com 4
povoações
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
76 / 164
Dados vectoriais em R: manipular SpatialPointsDataFrame
Caso se seleccione uma coluna, o resultado é também um objecto
SpatialPointsDataFrame cuja tabela de atributos tem apenas
essa coluna:
1
pov . e t r s [ , " designa " ] # devolve S p a t i a l P o i n t s D a t a F r a m e com um
ú n i c o a t r i b u t o ( designa )
O acesso às coordenadas dos pontos em
SpatialPointsDataFrame faz-se com a função coordinates:
1
c o o r d i n a t e s ( pov . e t r s )
Por exemplo, pode usar-se as coordenadas para representar a
localizações sobre uma imagem:
1
2
3
p l o t (mde . e t r s , x a x t = " n " , y a x t = " n " , box=FALSE , axes=FALSE , legend
=FALSE , c o l = t e r r a i n . c o l o r s ( 1 0 0 ) )
polygon ( l i m i t e . ap )
t e x t ( c o o r d i n a t e s ( pov . e t r s ) , as . c h a r a c t e r ( pov . etrs@data$ designa )
)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
77 / 164
Dados vectoriais em R: manipular SpatialPointsDataFrame
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
78 / 164
Dados vectoriais em R: criar SpatialPolygonsDataFrame
No exemplo seguinte, vamos transformar o polígono do limite da área
protegida num objecto SpatialPolygonsDataFrame.
A estrutura vai ter um único objecto multi-parte, que por sua vez tem
uma única parte, i.e. um único polígono. Primeiro, vai ser criada essa
parte a partir da matriz com duas colunas de coordenadas
limite.ap usando a função Polygon:
1
p a r t 1 <− Polygon ( l i m i t e . ap , h o l e =FALSE )
A opção hole=FALSE indica que é um polígono “positivo”, e não um
polígono a ser retirado a outra parte. Em seguida, cria-se o objecto
multi-parte combinando as partes, havendo neste caso apenas uma
parte, com a função Polygons:
1
mpart1 <− Polygons ( l i s t ( p a r t 1 ) , "SC" ) # c r i a r o b j e c t o m u l t i −
parte
A designação "SC" é arbitrária mas obrigatória: todos os objectos
multi-parte devem ter uma designação.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
79 / 164
Dados vectoriais em R: criar SpatialPolygonsDataFrame
O passo seguinte é criar a estrutura espacial completa do objecto
SpatialPolygonsDataFrame. Para tal combinam-se as
multi-partes e escolhe-se um CRS com a função SpatialPolygons:
1
ap . sp <− S p a t i a l P o l y g o n s ( l i s t ( mpart1 ) , p r o j 4 s t r i n g =CRS( e t r s ) )
Finalmente, associa-se ao objecto SpatialPolygons criado, a
tabela de atributos, representada por um objecto data.frame.
Os nomes das linhas de data.frame devem ser iguais aos nomes
dos objectos multi-parte que compõem a estrutura
SpatialPolygons. Neste caso, df tem uma única coluna nome e
uma única linha dado que existe apenas uma multi-parte:
1
d f <− data . frame ( nome=c ( " S i n t r a Cascais " ) , row . names=c ( "SC" ) )
Associar a componente geográfica com a tabela de atributos faz-se
com a função SpatialPolygonsDataFrame:
1
ap . e t r s <− SpatialPolygonsDataFrame ( ap . sp , d f )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
80 / 164
Dados vectoriais em R: estrutura de SpatialPolygonsDataFrame
As componentes de SpatialPolygonsDataFrame correspondem a
slots @data, @polygons, @plotOrder, @bbox, @proj4string. A
estrutura de SpatialPolygons é a seguinte:
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
81 / 164
Dados vectoriais em R: exportar objectos sp
Como foi visto anteriormente, pode exportar-se um cdg raster para
um ficheiro GeoTIFF com a função writeRaster:
1
w r i t e R a s t e r (mde . e t r s , f i l e = "mde_ s i n t r a _ c a s c a i s . t i f " , o v e r w r i t e =
TRUE)
Os objectos sp das classes SpatialPointsDataFrame,
SpatialLinesDataFrame e SpatialPolygonsDataFrame
podem ser exportados como um conjunto de ficheiros shapefile com a
função writeOGR:
1
2
writeOGR ( ap . e t r s , dsn=getwd ( ) , l a y e r = " L i m i t e " , d r i v e r = " ESRI
S h a p e f i l e " , o v e r w r i t e _ l a y e r =TRUE)
writeOGR ( pov . e t r s , dsn=getwd ( ) , l a y e r = " L o c a l i d a d e s " , d r i v e r = " ESRI
S h a p e f i l e " , o v e r w r i t e _ l a y e r =TRUE)
Sugestão: abrir am ArcGIS ou QGIS os ficheiros criados com os
comandos acima. Verificar os CRS dos conjuntos de dados e ver
ficheiro .prj para o caso dos dados vectoriais.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
82 / 164
Dados vectoriais em R: ler shapefile em R
Exercício com a totalidade dos dados de áreas protegidas: Ler
ficheiros do tipo shapefile para o conjunto de areas protegidas (fonte:
ICNF):
1
i c n f<−readOGR ( dsn=getwd ( ) , l a y e r = "AP_JUL_2014 " , encoding= " ISO8859
−1" )
O output icnf é o cdg SpatialPolygonsDataFrame. Qual é o
CRS desse cdg?
1
icnf@proj4string
Quantos multi-poligonos tem esse cdg?
1
l e n g t h ( icnf@polygons )
Como são as primeiras linhas da tabela de atributos?
1
head ( icnf@data )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
83 / 164
Dados vectoriais em R: selecção por atributos em objectos sp
Procurar a linha da tabela de atributos que tem Montejunto no atributo
Nome e verificar resultado;
1
2
i n d i c e<−which ( g r e p l ( p a t t e r n = " Montejunto " , icnf@data $NOME, i g n o r e .
case = TRUE) ) # í n d i c e é 12
icnf@data [ i n d i c e , ]
A função lógica grepl permite determinar se a sequência de
caractéres “Montejunto” pertence a alguns dos elementos do vector
icnf@data$NOME. O argumento ignore.case permite ignorar a
diferenças entre maiúsculas e minúsculas.
Para seleccionar essa 12o multi-parte de icnf selecciona-se a
respectiva linha tal como se faz para seleccionar uma linha da tabela
de atributos:
1
p l o t ( i c n f [ 1 2 , ] ) # i c n f [ 1 2 , ] ainda é SpatialPolygonsDataFrame
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
84 / 164
Dados vectoriais em R: manipular SpatialPolygonsDataFrame
Extrair informação de SpatialPolygonsDataFrame.
Matriz com a extensão de incf[12, ]:
1
i c n f [ 1 2 , ] @bbox # m a t r i x com l i n h a s " x " e " y " e colunas " min " e
"max"
ID da 12a multi-parte:
1
icnf@polygons [ [ 1 2 ] ] @ID
Área da 12a multi-parte:
1
icnf@polygons [ [ 1 2 ] ] @area
Centróide da 12a multi-parte:
1
icnf@polygons [ [ 1 2 ] ] @labpt
Valor do atributo "NOME" para a 12a multi-parte:
1
2
icnf@data [ 1 2 , "NOME" ] # ou , equivalentemente ,
i c n f [ 1 2 , ] @data$NOME
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
85 / 164
Dados vectoriais em R: manipular SpatialPolygonsDataFrame
Quantas partes tem a 12a multi-parte de icnf?
1
l e n g t h ( icnf@polygons [ [ 1 2 ] ] @Polygons ) #tem 6 p a r t e s
Qual é o tipo de cada parte: é uma parte “positiva” hole==FALSE ou
negativa hole==TRUE?
1
f o r ( i i n 1 : 6 ) p r i n t ( icnf@polygons [ [ 1 2 ] ] @Polygons [ [ i ] ] @hole )
Qual é a área de cada parte da 12a multi-parte?
1
f o r ( i i n 1 : 6 ) p r i n t ( icnf@polygons [ [ 1 2 ] ] @Polygons [ [ i ] ] @area )
Quais são as coordenadas dos pontos que delimitam a 3a parte de
icnf[12,] (i.e. da 12a multi-parte de icnf)?
1
p o l 3 <− icnf@polygons [ [ 1 2 ] ] @Polygons [ [ 3 ] ] @coords # m a t r i z com 2
colunas
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
86 / 164
Dados vectoriais em R: manipular SpatialPolygonsDataFrame
Exercício: Construir uma imagem da Paisagem Protegida da Serra de
Montejunto, indicando a área do polígono “positivo” e as áreas dos
polígonos “negativos”:
1
2
3
4
5
6
7
8
9
10
11
12
plot ( icnf [12 ,])
t e x t ( x= i c n f [ 1 2 , ] @bbox [ " x " , " min " ] ,
y= i c n f [ 1 2 , ] @bbox [ " y " , "max " ] ,
as . c h a r a c t e r ( icnf@data [ 1 2 , "NOME" ] ) , pos =4 , cex = . 9 )
f o r ( i i n 1 : l e n g t h ( icnf@polygons [ [ 1 2 ] ] @Polygons ) )
{
aux<−icnf@polygons [ [ 1 2 ] ] @Polygons [ [ i ] ] ;
i f ( aux@hole ) polygon ( aux@coords , c o l = " y e l l o w " )
t e x t ( x=aux@labpt [ 1 ] ,
y=aux@labpt [ 2 ] ,
paste ( round ( aux@area / 10000 ,1) , " ha " , sep= " " ) , cex = . 7 ,
pos =4)
}
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
87 / 164
Dados vectoriais em R: manipular SpatialPointsDataFrame
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
88 / 164
Dados vectoriais em R: manipular SpatialPolygonsDataFrame
Exercício: Escrever uma função para saber quantas partes de cada
tipo (“hole” ou não) há em cada multi-parte; o output é uma lista com
duas componentes: $positivos, vector de número de partes em
que hole=FALSE; $negativos, idem para hole=TRUE.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
t i p o . p a r t e s<− f u n c t i o n ( s p d f )
{
pos<−c ( ) ; neg<−c ( )
f o r ( i i n 1 : l e n g t h ( spdf@polygons ) ) # m u l t i −p a r t e s
{
p<−0 ; n<−0
f o r ( j i n 1 : l e n g t h ( spdf@polygons [ [ i ] ] @Polygons ) ) # p a r t e s
{
i f ( spdf@polygons [ [ i ] ] @Polygons [ [ j ] ] @hole ) n<−n+1 e l s e p
<−p+1 # incrementa n se n e g a t i v a ; senão incrementa p
}
pos<−c ( pos , p ) ; neg<−c ( neg , n )
}
r e t u r n ( l i s t ( p o s i t i v o s =pos , n e g a t i v o s =neg ) )
}
t i p o . p a r t e s ( i c n f ) # devolve r e s u l t a d o para dados do ICNF
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
89 / 164
Dados vectoriais em R: manipular SpatialPolygonsDataFrame
Exercício mais complexo: Um investigador entende que é importante
definir a zona de influência sobre cada área protegida como uma
elipse que envolve a área protegida, sem ter em conta as partes
“negativas”. Pretende-se então construir um
SpatialPolygonsDataFrame designado por elipses que
contenha essa informação.
Cada elipse pode ser criada com o package ellipse que devolve
uma matriz com colunas “x” e “y” em torno de um centro (argumento
center) e com dimensões definidas por uma matriz, que neste caso
é a matriz de variância-covariância do limite de cada polígono.
1
library ( ellipse )
O objecto multipartes é uma lista que vai armazenar as várias
elipses:
1
m u l t i p a r t e s<− l i s t ( )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
90 / 164
Dados vectoriais em R: manipular SpatialPolygonsDataFrame
O ciclo principal passa por todas as partes, e para as que sejam
“positivas” constroi um objecto Polygon usando o polígono devolvido
pela função ellipse. Nota: a função round é apenas para
assegurar que o primeiro e o último ponto da elipse coincidem.
1
2
3
4
5
6
7
8
9
10
11
12
13
f o r ( i i n 1 : l e n g t h ( icnf@polygons ) )
{
p a r t e s<− l i s t ( ) ; np<−0
f o r ( j i n 1 : l e n g t h ( icnf@polygons [ [ i ] ] @Polygons ) )
{
i f ( icnf@polygons [ [ i ] ] @Polygons [ [ j ] ] @hole==FALSE )
{
np<−np+1
p a r t e s [ [ np ] ] <− Polygon ( round ( e l l i p s e ( cov ( icnf@polygons [ [
i ] ] @Polygons [ [ j ] ] @coords ) , c e n t r e = icnf@polygons [ [ i ] ]
@Polygons [ [ j ] ] @labpt ) , 2 ) , h o l e =FALSE )
}
}
m u l t i p a r t e s [ [ i ] ] <−Polygons ( p a r t e s , ID=icnf@polygons [ [ i ] ] @ID )
}
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
91 / 164
Dados vectoriais em R: manipular SpatialPolygonsDataFrame
Para concluir:
1
Define-se o CRS, igual ao de icnf;
2
Dado que os atributos de elipses são os mesmos do que os de
icnf, usa-se simplesmente a tabela de atributos de icnf para
definir o objecto SpatialPolygonsDataFrame:
Nota: no passo anterior usou-se como ID de cada multi-parte os IDs
de icnf. Assim, os IDs e os nomes das linhas da tabela de atributos
são iguais como exigido por SpatialPolygonsDataFrame.
1
2
e l i p s e s . e t r s<−S p a t i a l P o l y g o n s ( m u l t i p a r t e s , p r o j 4 s t r i n g =
icnf@proj4string )
e l i p s e s <− SpatialPolygonsDataFrame ( e l i p s e s . e t r s , data=icnf@data
)
Finalmente, o resultado pode ser exportado como uma shapefile:
1
writeOGR ( e l i p s e s , dsn=getwd ( ) , l a y e r = " E l i p s e s " , d r i v e r = " ESRI
S h a p e f i l e " , o v e r w r i t e _ l a y e r =TRUE)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
92 / 164
Dados vectoriais em R: representar elipses e áreas protegidas
sobre dados de concelhos
A shapefile designada conc_1998, obtida no site do Instituto do
Ambiente, contém a representação dos concelhos de Portugal
Continental em 1998. Pretende-se sobrepor as áreas protegidas e as
elipses sobre esse cdg.
1
2
3
4
5
concelhos<−readOGR ( dsn=getwd ( ) , l a y e r = " conc_1998 " , encoding= "
ISO8859−1" ) # cuidado com o CRS
c o n c e l h o s @ p r o j 4 s t r i n g # v e r i f i c a r que não tem t ra n sf o rm a çã o de
datum ! !
igeoe . g r i d <− " + p r o j =tmerc + l a t _0=39.66666666666666 + n a d g r i d s =
ptLX_e89 . gsb + l o n _0=1 +k=1 +x_0=200000 +y_0=300000 + e l l p s =
i n t l +pm= l i s b o n + u n i t s =m"
c o n c e l h o s @ p r o j 4 s t r i n g<−CRS( igeoe . g r i d )
concelhos . e t r s<−spTransform ( concelhos ,CRS( e t r s ) )
Nota: para usar o método das grelhas, o ficheiro ptLX_e89.gsb deve
estar na pasta dada por
system.file("proj", package = "rgdal").
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
93 / 164
1
2
3
Dados vectoriais em R: representar elipses e áreas protegidas
sobre dados de concelhos
Agora, todos os cdg têm o mesmo
sistema de coordenadas e podem
ser sobrepostos:
p l o t ( concelhos . e t r s )
p l o t ( e l i p s e s , c o l = " green " , add=
TRUE)
p l o t ( i c n f , c o l = " red " , add=TRUE)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
94 / 164
Extrair dados matriciais usando dados vectoriais como
delimitadores
Para ilustrar este procedimento vai usar-se um conjunto de parcelas
agrícolas no Ribatejo para as quais é conhecida a cultura. Desse cdg
foram extraídas as parcelas com ocupação de sorgo e com ocupação
de luzerna: o cdg resultante está disponível na shapefile
SorgoLuzerna:
1
2
3
sorgo . l u z e r n a<− readOGR ( dsn=getwd ( ) , l a y e r = " SorgoLuzerna " ,
encoding= " ISO8859−1" )
# C r i a r v e c t o r com os nomes das c u l t u r a s para cada p a r c e l a
nomes . p a r c e l a s<−as . c h a r a c t e r ( sorgo . luzerna@data$ c u l t u r a )
Uma imagem multiespectral Landsat pode ser acedida com:
1
2
3
f i c h s <− l i s t . f i l e s ( p a t t e r n = " banda " )
s <− s t a c k ( as . l i s t ( f i c h s ) )
b <− b r i c k ( s )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
95 / 164
1
2
3
4
5
6
Extrair dados matriciais usando dados vectoriais como
delimitadores
Como os cdg têm o mesmo
sistema de coordenadas, podem
ser sobrepostos. Aqui usa-se a
extensão do cdg vectorial para
delimitar plotRGB:
# v e r i f i c a r CRS
sorgo . l u z e r n a @ p r o j 4 s t r i n g
b@crs
# c o n s t r u i r imagem
plotRGB ( b , r =4 ,g=3 ,b=2 , s t r e t c h = "
l i n " , e x t =sorgo . luzerna@bbox
)
p l o t ( sorgo . l u z e r n a , add=TRUE, c o l
=" yellow " )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
96 / 164
Extrair dados matriciais usando dados vectoriais
Para fazer a extração das células do objecto RasterBrick usa-se
função extract que devolve os valores das células no interior de
cada polígono do objecto sorgo.luzerna. Neste exemplo, o output
é uma tabela (data.frame). O argumento weights permite que seja
devolvida a proporção de cada célula de b que está contida em cada
polígono de sorgo.luzerna. Nota: quando weights=TRUE o
tempo de execução é muito mais elevado.
Aplicar extract usando a+enas a primeira parcela de
sorgo.luzerna para compreender melhor o resultado:
1
2
o u t<− e x t r a c t ( b , sorgo . l u z e r n a [ 1 , ] , w e i g h t s =TRUE, normalizeWeights =
FALSE , c e l l n u m b e r s =TRUE)
outmat<−o u t [ [ 1 ] ]
Construir uma imagem que indica a proporção de cada célula (pixel)
que pertence à parcela:
1
2
3
4
plotRGB ( b , r =4 ,g=3 ,b=2 , s t r e t c h = " l i n " , e x t =sorgo . l u z e r n a [ 1 , ] @bbox )
p l o t ( sorgo . l u z e r n a [ 1 , ] , add=TRUE, c o l = " y e l l o w " )
xy<−c o o r d i n a t e s ( b ) [ outmat [ , " c e l l " ] , ]
tManuel
e x t (Campagnolo
xy [ , 1 ] (ISA/ULisboa)
, xy [ , 2 ] , outmat [ , " wSIGs
e i gcom
h t R" ] )
5 a 10 de Fevereiro de 2015
97 / 164
Extrair dados matriciais usando dados vectoriais
A extracção para a totalidade das parcelas devolve uma lista; cada
elementos da lista é uma matriz como a matriz outmat do exemplo
acima.
1
2
o u t<− e x t r a c t ( b , sorgo . l u z e r n a , w e i g h t s =TRUE, normalizeWeights =
FALSE )
length ( out )
Pode “empilhar-se” todas as matrizes para obter uma única matriz
com todos os valores extraídos em todas as parcelas fazendo o
seguinte ciclo:
1
2
3
outmat<−o u t [ [ 1 ] ]
f o r ( i i n 2 : l e n g t h ( o u t ) ) outmat<−r b i n d ( outmat , o u t [ [ i ] ] )
dim ( outmat ) # devolve o tamanho da m a t r i z r e s u l t a n t e
Para seleccionar as células que estão contidas a 95% nas parcelas:
1
o u t c e l <− outmat [ outmat [ , " w e i g h t " ] > 0 . 9 5 , ]
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
98 / 164
Conjuntos de dados geográficos vectoriais em
R: análise espacial com o package rgeos
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
99 / 164
Dados vectoriais: conceitos fundamentais de análise espacial
O OGC OpenGIS Implementation Standard for Geographic information
/ ISO 19125 define métodos para objectos geométricos:
Objectos geométricos do tipo point, line, polygon, multi-point, etc,
associados a um sistema de coordenadas (ver Slide 72);
Métodos sobre objectos geométricos devolvem propriedades tais
como dimensão, fronteira, área, centroide, etc (ver
Slide 101 );
Métodos para testar relações espaciais entre objectos geométricos
equals, disjoint, intersects, touches, crosses, within,
contains, overlaps and relate, com valor lógico TRUE ou
FALSE (ver Slide 106);
Métodos que suportam análise espacial distância, que devolve
uma distância, e buffer, convex hull, intersection, union,
difference, symmetric difference, que devolvem novos
objectos geométricos (ver Slide 114).
Fonte: www.opengeospatial.org/standards/sfa
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
100 / 164
Métodos sobre objectos geométricos
Objectos do tipo “ponto” (p) têm dimensão 0,
Objectos do tipo “linha” (L) têm dimensão 1,
Objectos do tipo “polígono” (P) têm dimensão 2
Os objectos geométricos têm Interior, Fronteira e Exterior:
1
Se p é do tipo ponto, então I(p) = p, F(p) é vazio, e E(p) são
todos os pontos que não estão em p;
2
Se L é do tipo linha, então F(L) são os extremos da linha, I(L)
são todos os pontos da linha excepto os extremos, e o exterior
E(L) são os todos os outros pontos;
3
Se P é do tipo polígono, F(P)é a fronteira, I(P) é o conjunto de
pontos de P que não estão na fronteira e E(P) são os restantes
pontos.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
101 / 164
Some basic properties of geometric objects: area and centroid
Geometric objects of dimension 2 have an area, which can be
determined with an appropriate function.
The area is measured with respect to the Coordinate Reference
System of the object and therefore will depend on the map projection
and on the distortions associated to the location.
Among other basic properties, 2-dimensional geometric objects also
have a centroid, which is the point of the center of mass of the object.
However, the centroid might be located outside the geometric object.
Hence, the OGC standard also defines the method PointOnSurface
which returns a point that lies on the object.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
102 / 164
Análise espacial sobre objectos geométricos: o package rgeos
O package rgeos contém funções para implementar todos os
métodos da norma descrita atrás. Para exemplificar vamos usar dados
analisados anteriormente (dados do ICNF):
1
2
l i b r a r y ( rgeos )
i c n f<−readOGR ( dsn=getwd ( ) , l a y e r = "AP_JUL_2014 " , encoding= " ISO8859
−1" )
Para usar funções do package rgeos, as geometrias dos objectos
devem ser válidas; a função gIsvalid do package rgeos permite
saber quais são os objectos geométricos válidos:
1
2
v a l i d<− g I s V a l i d ( i c n f , b y i d =TRUE, reason=TRUE)
which ( v a l i d ! = " V a l i d Geometry " ) # 13 55 58 61 i n v á l i d o s
Vamos restringir a análise aos objectos geométricos válidos:
1
i c n f v<− i c n f [ which ( v a l i d == " V a l i d Geometry " ) , ]
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
103 / 164
Análise espacial: métodos para objectos geométricos
Determinar áreas as áreas protegidas com a função gArea:
1
2
gArea ( i c n f v ) # devolve a área t o t a l
gArea ( i c n f v , b y i d =TRUE) # áreas de todos os p o l í g o n o s
A função gPointOnSurface aplicada a um objecto espacial sp
devolve um objecto SpatialPoints. Cada elemento do objecto de
input é intersectado por um ponto do output. Esta função é útil por
exemplo para colocar texto sobre a imagem do input. Nota: o slot
@labpt visto anteriormente devolve o centróide, que pode não
intersectar o objecto correspondente.
1
2
3
4
p l o t ( i c n f v , c o l = " green " , b o r d e r =FALSE )
areas <− round ( as . v e c t o r ( gArea ( i c n f v , b y i d =TRUE) ) / 10000) # areas
em ha
xy <− c o o r d i n a t e s ( gPointOnSurface ( i c n f v , b y i d =TRUE) )
t e x t ( xy [ areas > 2 0 0 0 0 , 1 ] , xy [ areas >20000 ,2] , paste ( areas [ areas
>20000]) , cex = 0 . 7 )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
104 / 164
Análise espacial: métodos para objectos geométricos
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
105 / 164
Métodos para testar relações espaciais entre objectos
geométricos
A norma OGC/ISO (ver Slide 100) define um conjunto de operadores
com valor lógico (TRUE ou FALSE ) para testar relações espaciais
entre objectos geométricos ponto, Linha, Polígono.
1
2
3
4
5
6
7
8
9
Equal: os objectos coincidem;
Disjoint: os objectos não têm pontos em comum;
Touches aplica-se a 2 objectos desde que algum tenha dimensão
maior que 0: a Touches b se as fronteiras se intersectam;
Para a/b de tipo p/L, p/P, L/L or L/P, a Crosses b se os objectos
se intersectam mas nenhum contém o outro;
a Within b se a não intersecta o exterior de b;
Overlaps aplica-se a objectos com mesma dimensão: é TRUE se
se os objectos se intersectam mas nenhum contém o outro;
a Contains b se b Within a;
a Intersects b se a e b têm algum ponto em comum;
Relates permite definir qualquer relação espacial.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
106 / 164
Relações espaciais entre objectos geométricos (ilustração)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
107 / 164
Métodos para testar relações espaciais entre objectos
geométricos com package rgeos
Estes testes permitem selecções por localização. A ideia é selecionar
um subconjunto de objectos de uma “feature class” com base na
relação espacial relativamente a outra “feature class”.
No exemplo abaixo seleccionam-se os multi-polígonos (concelhos)
que intersectam as linhas do cdg rios que têm “Rio Mondego” na
designação:
1
2
3
4
5
6
7
nomerio<−" r i o mondego "
r i o<− r i o s . e t r s [ which ( g r e p l ( p a t t e r n =nomerio , r i o s . etrs@data$
DESIGNACAO, i g n o r e . case = TRUE) ) , ]
# i t s é uma m a t r i z l ó g i c a (TRUE/ FALSE ) 49 * 282
# há 49 m u l t i −l i n h a s que compõem r i o
# há 282 m u l t i −p o l í g o n o s que compõem concelhos . e t r s
i t s <− g I n t e r s e c t s ( concelhos . e t r s , r i o , b y i d =TRUE)
dim ( i t s )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
108 / 164
Métodos para testar relações espaciais entre objectos
geométricos com package rgeos
Para selecionar os concelhos que são intersectados pelo rio é
necessário saber quais as colunas da matriz its que têm pelo menos
um TRUE:
1
2
# d e f i n i r v e c t o r de TRUE e FALSE para os concelhos
aux<−as . l o g i c a l ( a p p l y ( i t s , 2 , max ) )
Os concelhos seleccionados formam o objecto
concelhos.etrs[aux,]. Por exemplo, podemos criar um mapa
com esses concelhos:
1
2
3
p l o t ( concelhos . e t r s , x l i m =c ( −95000 , 151000) , y l i m =c
(11000 ,136000) )
p l o t ( concelhos . e t r s [ aux , ] , add=TRUE, c o l = " y e l l o w " )
p l o t ( r i o , add=TRUE, c o l = " b l u e " )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
109 / 164
Métodos para testar relações espaciais entre objectos
geométricos com package rgeos
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
110 / 164
Métodos para testar relações espaciais entre objectos
geométricos com package rgeos
Vamos agora identificar os concelhos no Parque da Serra da Estrela.
Em primeiro lugar vamos seleccionar (selecção por atributos) esse
parque:
1
2
nomeap <−" e s t r e l a "
parque<− i c n f v [ which ( g r e p l ( p a t t e r n =nomeap , icnfv@data $NOME, i g n o r e
. case = TRUE) ) , ]
Seleccionar concelhos que intersectam o parque:
1
2
3
i t s <− g I n t e r s e c t s ( concelhos . e t r s , parque , b y i d =TRUE)
aux<−as . l o g i c a l ( a p p l y ( i t s , 2 , max ) ) # v e c t o r l o g i c o
concelhos . e t r s [ aux , ]
Seleccionar concelhos que que estão contidos no parque:
1
2
3
i t s <−g W i t h i n ( concelhos . e t r s , parque , b y i d =TRUE)
aux<−as . l o g i c a l ( a p p l y ( i t s , 2 , max ) ) # v e c t o r l o g i c o
concelhos . e t r s [ aux , ]
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
111 / 164
Métodos para testar relações espaciais entre objectos
geométricos com package rgeos
Para criar mapa:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
p l o t ( parque , c o l = " green " )
# concelhos que i n t e r s e c t a m o Parque
i t s <− g I n t e r s e c t s ( concelhos . e t r s , parque , b y i d =TRUE)
aux<−as . l o g i c a l ( a p p l y ( i t s , 2 , max ) ) # v e c t o r l o g i c o
p l o t ( concelhos . e t r s [ aux , ] , add=TRUE, b o r d e r = " red " )
# e s c r e v e r nomes dos concelhos s e l e c c i o n a d o s
xy<−c o o r d i n a t e s ( gPointOnSurface ( concelhos . e t r s [ aux , ] , b y i d =TRUE)
)
t e x t ( xy [ , 1 ] , xy [ , 2 ] , concelhos . etrs@data [ aux , "CONCELHO" ] , cex = 0 . 7 )
# concelhos t o t a l m e n t e c o n t i d o s no parque :
i t s <−g W i t h i n ( concelhos . e t r s , parque , b y i d =TRUE)
aux<−as . l o g i c a l ( a p p l y ( i t s , 2 , max ) ) # v e c t o r l o g i c o
p l o t ( concelhos . e t r s [ aux , ] , add=TRUE, c o l = " red " )
# e s c r e v e r nome do concelho s e l e c c i o n a d o
xy<−c o o r d i n a t e s ( gPointOnSurface ( concelhos . e t r s [ aux , ] , b y i d =TRUE)
)
t e x t ( xy [ , 1 ] , xy [ , 2 ] , concelhos . etrs@data [ aux , "CONCELHO" ] , cex = 0 . 7 )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
112 / 164
Métodos para testar relações espaciais entre objectos
geométricos com package rgeos
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
113 / 164
Métodos que suportam análise espacial
A norma OGC/ISO (ver Slide 100) define operadores que devolvem a
distância entre objectos geométricos ou novos objectos geométricos
de acordo com alguma relação espacial. Dados objectos a e b:
1
Distance devolve a distância mais curta entre pontos de a e de b;
2
Buffer de a coom distância d devolve o objecto geométrico cujos
pontos estão uma distância menor ou igual a d de a;
3
ConvexHull de a devolve um objecto geométrico que é o invólucro
convexo de a;
4
Intersection devolve o objecto geométrico cujos pontos estão em
a e em b;
5
Union devolve o objecto geométrico cujos pontos estão em a ou
em b;
6
Difference entre a e b devolve o objecto geométrico com os
pontos que estão em a mas não em b (operação não comutativa);
7
Symmetric Difference entre a e b devolve o objecto geométrico
cujos pontos satisfazem apenas uma das condições: (1) estão
em a mas não em b; (2) estão em b mas não em a.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
114 / 164
Métodos que suportam análise espacial com package rgeos
A função gDistance permite calcular distências entre “features” de
dois objectos sp. No exercício seguinte determinam-se as distâncias
entre os concelhos e as áreas protegidas com área superior a 1000ha.
O resultado da função é uma matriz de distâncias.
Em primeiro lugar define-se um novo cdg com as áreas protegidas
com área superior a 1000ha:
1
2
areas <− round ( as . v e c t o r ( gArea ( i c n f v , b y i d =TRUE) ) / 10000)
i c n f g<− i c n f v [ areas >1000 ,]
Em seguida calculam-se as distâncias; d é uma matriz 27 × 282 que
contem as distâncias entre concelhos e APs:
1
2
3
d<−gDistance ( concelhos . e t r s , i c n f g , b y i d =TRUE) # demora uns 2mn a
calcular
colnames ( d )<−as . c h a r a c t e r ( concelhos . etrs@data$CONCELHO)
rownames ( d )<−as . c h a r a c t e r ( icnfg@data $NOME)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
115 / 164
Métodos que suportam análise espacial com package rgeos
Qual é distância mínima (em Km) de cada concelho a uma área
protegida com mais de 1000ha?
1
dc<−round ( a p p l y ( d , 2 , min ) / 1000)
Qual o concelho que fica à maior distância?
1
dc [ dc==max ( dc ) ] # Mira , a 75 km
Qual é a área protegida mais próxima desse concelho
1
2
3
i n d<−which ( dc==max ( dc ) ) # devolve o í n d i c e do concelho ( coluna
de d )
aux<−d [ , i n d ] # dá as d i s t â n c i a s desse concelho a todas as áreas
protegidas
aux [ aux==min ( aux ) ] # devolve a área p r o t e g i d a à menor d i s t â n c i a
do concelho
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
116 / 164
1
2
3
Métodos que suportam análise espacial com package rgeos
Fazer figura:
p l o t ( concelhos . e t r s )
p l o t ( concelhos . e t r s [ ind , ] , c o l = "
red " , add=TRUE)
p l o t ( i c n f g , c o l = " green " , add=TRUE
)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
117 / 164
Métodos que suportam análise espacial com package rgeos
Uma operação fundamental em SIGs é a definição de buffers. Em
rgeos a funçao respectiva é gBuffer.
Construir um buffer em torno da área protegida da Serra de
Montejunto:
1
2
p l o t ( g B u f f e r ( i c n f v [ 1 2 , ] , b y i d =TRUE, w i d t h =100) , b o r d e r = " red " )
p l o t ( i c n f v [ 1 2 , ] , add=TRUE)
Nota: a função gBuffer aceita argumentos adicionais que são úteis
sobretudo para buffers de linhas:
1
quadsegs: Number of line segments to use to approximate a
quarter circle.
2
capStyle: Style of cap to use at the ends of the geometry. Allowed
values: ROUND,FLAT,SQUARE
3
joinStyle: Style to use for joints in the geometry. Allowed values:
ROUND,MITRE,BEVEL
4
mitreLimit: Numerical value that specifies how far a joint can
extend if a mitre join style is used.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
118 / 164
Métodos que suportam análise espacial com package rgeos
Pode, obviamente, criar-se vários anéis de buffer com um ciclo. Neste
caso usam-se cores do package abaixo para os anéis.
1
2
l i b r a r y ( RColorBrewer ) # para c r i a r p a l e t a s
d i s p l a y . brewer . a l l ( ) # v e r p o s s í v e i s p a l e t a s de cores
O ciclo pode ser o seguinte:
1
2
3
p l o t ( g B u f f e r ( i c n f v [ 1 2 , ] , w i d t h =1100) , b o r d e r =FALSE )
f o r ( i i n 1 0 : 1 ) p l o t ( g B u f f e r ( i c n f v [ 1 2 , ] , b y i d =TRUE, w i d t h =100 * i ) ,
add=TRUE, c o l =brewer . p a l ( n=12 , " S p e c t r a l " ) [ i ] )
p l o t ( i c n f v [ 1 2 , ] , add=TRUE, c o l = " green " )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
119 / 164
Métodos que suportam análise espacial com package rgeos
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
120 / 164
Métodos que suportam análise espacial com package rgeos
A dissolução de “features” faz-se com a função gUnaryUnion.
Para exemplificar, criar SpatialPolygonsDataFrame dos distritos de
Portugal Continental:
1
2
3
4
# 1o : d e f i n i r um v e c t o r de IDs para agrupar as m u l t i −p a r t e s do
input
i d . d i s t r i t o s <−as . c h a r a c t e r ( concelhos . etrs@data [ , " DISTRITO " ] )
# 2o : usar esse v e c t o r como argumento de gUnaryUnion
d i s t r i t o s <−gUnaryUnion ( concelhos . e t r s , i d = i d . d i s t r i t o s )
gUnaryUnion devolve um objecto SpatialPolygons sem tabela de
atributos. Para obter um SpatialPolygonsDataFrame temos que criar a
tabela de atributos.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
121 / 164
Métodos que suportam análise espacial com package rgeos
Construir a tabela de atributos e associá-la ao objecto
SpatialPolygons:
1
2
3
4
# IDs associados às m u l t i −p a r t e s de d i s t r i t o s
i d s<−c ( ) ; f o r ( i i n 1 : 1 8 ) i d s<−c ( i d s , d i s t r i t o s @ p o l y g o n s [ [ i ] ] @ID
)
# ou , bem mais simples ,
names ( d i s t r i t o s )
Construir um vector de dimensão 18 com a região respectiva:
1
2
3
4
5
6
7
d f <− concelhos . etrs@data
r e g i o e s<−c ( ) ;
f o r ( i in 1:18)
{
r g s<−as . c h a r a c t e r ( d f [ as . c h a r a c t e r ( d f $DISTRITO ) == i d s [ i ] , "
REGIAO" ] )
r e g i o e s [ i ]<−r g s [ 1 ]
}
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
122 / 164
1
2
3
Métodos que suportam análise espacial com package rgeos
Criar data.frame para associar a
distritos:
d f . d i s t r i t o s <−data . frame ( nome= i d s ,
r e g i a o = r e g i o e s , row . names= i d s )
d i s t r i t o s . e t r s<−
SpatialPolygonsDataFrame (
d i s t r i t o s , data= d f . d i s t r i t o s )
p l o t ( d i s t r i t o s . e t r s , c o l =brewer . p a l ( n
=10 , " S p e c t r a l " ) )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
123 / 164
Métodos que suportam análise espacial com package rgeos
A operação espacial de intersecção de "features" de objectos sp
realiza-se com a função gIntersection. Voltando ao exemplo dos
concelhos no Parque Natural da Serra da Estrela:
1
o u t<− g I n t e r s e c t i o n ( parque , concelhos . e t r s , b y i d =TRUE)
De novo, out é de classe SpatialPolygons. Para obter um objecto
de classe SpatialPolygonsDataFrame é preciso definir a tabela
de atributos. Pretende-se ter na tabela de atributos o nome do
concelho, do distrito, e o nome do Parque.
Os IDs associados às multi-partes de out são:
1
names ( o u t ) # "34 115" "34 119" "34 120" "34 126" "34 131" "34
136" "34 139" "34 143"
Os IDs dos multi-partes são definidos pela concatenação dos IDs de
parque, neste caso sempre 34, e dos IDs de concelhos.etrs.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
124 / 164
Métodos que suportam análise espacial com package rgeos
Para saber qual é o parque e qual é o concelho é preciso
1
extrair o ID respectivo;
2
usar as tabelas de atributos de parque e concelhos.etrs para
obter as designações.
A função strsplit separa e devolve uma lista com 8 componentes:
cada componente da lista é um vector com dois elementos:
1
2
3
# o 1o elemento é r e l a t i v o a parque
# o 2o elemento é r e l a t i c o a concelhos . e t r s
s p l<− s t r s p l i t ( names ( o u t ) , s p l i t = " " )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
125 / 164
Métodos que suportam análise espacial com package rgeos
Criar tabela de atributos:
1
2
3
4
5
6
7
8
9
10
nomes . concelhos<−c ( )
nomes . d i s t r i t o s <−c ( )
nomes . ap<−c ( )
f o r ( i i n 1: length ( out ) )
{
nomes . ap<−c ( nomes . ap , as . c h a r a c t e r ( parque@data [ s p l [ [ i ] ] [ 1 ] , "
NOME" ] ) )
nomes . concelhos<−c ( nomes . concelhos , as . c h a r a c t e r ( concelhos .
etrs@data [ s p l [ [ i ] ] [ 2 ] , "CONCELHO" ] ) )
nomes . d i s t r i t o s <−c ( nomes . d i s t r i t o s , as . c h a r a c t e r ( concelhos .
etrs@data [ s p l [ [ i ] ] [ 2 ] , " DISTRITO " ] ) )
}
d f . o u t <− data . frame ( ap=nomes . ap , concelho=nomes . concelhos ,
d i s t r i t o =nomes . d i s t r i t o s , row . names=names ( o u t ) )
e associar essa tabela a out:
1
o u t . e t r s<−SpatialPolygonsDataFrame ( out , data= d f . o u t )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
126 / 164
1
2
3
4
5
6
7
Métodos que suportam análise espacial com package rgeos
Criar imagem:
p l o t ( o u t . e t r s , c o l =brewer . p a l ( n=
length ( out . e t r s ) , " Spectral
"))
f o r ( i i n 1: length ( out . e t r s ) )
{
xy<−o u t . etrs@polygons [ [ i ] ]
@labpt
t e x t ( xy [ 1 ] , xy [ 2 ] , o u t .
etrs@data [ i , " concelho " ] ,
pos =3 , cex = 0 . 5 )
t e x t ( xy [ 1 ] , xy [ 2 ] , o u t .
etrs@data [ i , " d i s t r i t o " ] ,
pos =1 , cex = 0 . 5 )
}
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
127 / 164
Métodos que suportam análise espacial com package rgeos
A função gConvexHull devolve o invólucro convexo:
1
2
p l o t ( i c n f v , c o l = " green " , b o r d e r =FALSE )
p l o t ( gConvexHull ( i c n f v , b y i d =TRUE) , add=TRUE)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
128 / 164
Dados matriciais: leitura e escrita com rgdal:
1
2
3
l i b r a r y ( rgdal )
# ler ficheiro
g <− readGDAL ( " landsat8pan . t i f " )
g é um objecto da classe SpatialGridDataFrame (package sp).
1
2
summary ( g ) # descreve g
s t r ( g ) # devolve a e s t r u t u r a de g
Os objectos das classes de sp têm os seguintes "slots":
1
@bbox: extensão geográfica dos dados;
2
@proj4string: sistema de coordenadas, de classe CRS;
3
@data: contém os valores associados às observações
geográficas.
Para a classe SpatialGridDataFrame, há o "slot" adicional:
1
@grid: estrutura espacial dos dados, i.e. o tamanho (resolução)
e coordenadas das células (ou pixels) do conjunto de dados
matricial.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
129 / 164
Dados matriciais: leitura e escrita com rgdal:
1
2
3
4
5
6
7
8
# g@data é uma das " s l o t s " do o b j e c t o : é uma data . frame
# os nomes das colunas de g@data são então dados por
names ( g@data ) # band1
# histograma dos v a l o r e s da imagem :
h i s t ( g@data$band1 )
cores <− gray ( seq ( 0 , 1 , l e n g t h . o u t =100) )
# image é melhor do que p l o t , que b l o q u e i a para l a n d s a t 8 pan
(225 * 10^6 p i x e l s )
image ( g , c o l =cores , z l i m = q u a n t i l e ( g@data$band1 , probs=c ( 0 , . 9 9 ) ) )
Exportar o objecto g para um ficheiro GeoTIFF:
1
writeGDAL ( g , fname= " o u t . t i f " ) #
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
130 / 164
Em jeito de conclusão
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
131 / 164
Interpolação espacial com R
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
132 / 164
Interpolação espacial
Definition (Spatial interpolation)
It is the interpolation of measurements at point locations to the whole
space. More precisely, given a set of measurements {z1 , . . . , zn } at
spatial locations {(x1 , y1 ), . . . , (xn , yn )}, spatial interpolation
determines a function z = f (x, y ) that interpolates the zi
measurements to all locations within the range of (xi , yi ) locations.
The point samples
can be elevation, temperature, precipitation
or other continuous
variables.
Figure: www.innovativegis.com/
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
133 / 164
Interpolação espacial em imagens por re-amostragem
Os dados de altimetria SRTM podem ser re-amostrados sobre uma
grelha de células com um maior número de células do que o original.
No exemplo que se segue vão ser usados dados de altimetria da
Serra da Estrela.
1
2
3
4
5
6
7
8
t i l e <−"N40W008"
# d e s c a r r e g a r dados de a l t i m e t r i a SRTM dessa r e g i ã o
URL<−paste ( " h t t p : / / dds . c r . usgs . gov / srtm / v e r s i o n 2 _1 /SRTM3/
E u r a s i a / " , t i l e , " . h g t . z i p " , sep= " " ) # SE
download . f i l e ( u r l =URL, d e s t f i l e = paste ( t i l e , " . h g t . z i p " , sep= " " ) ,
mode= "wb" ) # f a z download para pasta de t r a b a l h o
u n z i p ( z i p f i l e = paste ( t i l e , " . h g t . z i p " , sep= " " ) ) # f i c h e i r o . h g t
srtm<− r a s t e r ( paste ( t i l e , " . h g t " , sep= " " ) )
# r e c o r t a r na r e g i ã o das Penhas Douradas
srtm . crop<−crop ( srtm , e x t e n t ( c ( −7.62 , −7.54 ,40.32 ,40.4) ) )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
134 / 164
Interpolação espacial em imagens por re-amostragem
1
2
3
4
5
# c o n v e r t e r MDE para coordenadas ETRS
e t r s<−" + p r o j =tmerc + l a t _0=39.6682583 + l o n _0=−8.1331083 +k=1 +x_
0=0 +y_0=0 + e l l p s =GRS80 + u n i t s =m"
mde<−p r o j e c t R a s t e r ( srtm . crop , c r s = e t r s , method= " b i l i n e a r " )
# o b s e r v a r o MDE
p l o t ( mde , c o l =topo . c o l o r s ( 1 0 0 ) )
Realizar uma re-amostragem bi-linear para uma resolução mais fina
(20m):
1
2
3
4
5
6
7
mde20<−p r o j e c t R a s t e r ( srtm . crop , r e s =20 , c r s = e t r s , method= " b i l i n e a r
")
# nota : as extensões são d i f e r e n t e s :
mde@extent ;
mde20@extent
# para f i c a r exactamente com a mesma extensão
v a z i o<− r a s t e r ( e x t =mde@extent , c r s =mde@crs , r e s o l u t i o n =20)
mde20<−p r o j e c t R a s t e r ( from=srtm . crop , t o = vazio , method= " b i l i n e a r " )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
135 / 164
Interpolação espacial em imagens por re-amostragem
Adicionar à imagem do MDE curvas de nível e comparar resultados.
Para os dados de altimetria originais simplesmente re-projectados em
coordendas cartográficas:
1
2
3
4
l o c s<−as . v e c t o r ( mde@extent )
p l o t ( mde , c o l =topo . c o l o r s ( 1 0 0 ) , box=FALSE , axes=FALSE , legend=FALSE
)
c o n t o u r ( mde , add=TRUE)
t e x t ( x=mean ( l o c s [ 1 : 2 ] ) , y= q u a n t i l e ( l o c s [ 3 : 4 ] , . 9 ) , paste ( "
r e s o l u ç ã o = " , paste ( round ( r e s (mde) ) , c o l l a p s e = " , " ) ) , cex = 0 . 7 )
Para o cdg re-amostrado para uma resolução de 20 × 20 metros:
1
2
3
p l o t ( mde20 , c o l =topo . c o l o r s ( 1 0 0 ) , box=FALSE , axes=FALSE , legend=
FALSE )
c o n t o u r ( mde , add=TRUE)
t e x t ( x=mean ( l o c s [ 1 : 2 ] ) , y= q u a n t i l e ( l o c s [ 3 : 4 ] , . 9 ) , paste ( "
r e s o l u ç ã o = " , paste ( round ( r e s ( mde20 ) ) , c o l l a p s e = " , " ) ) , cex = 0 . 7 )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
136 / 164
Interpolação espacial em imagens por re-amostragem
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
137 / 164
Spatial interpolation methods
Different interpolation methods yield data of varying quality and can
strongly influence modeling results.
linear, e.g. interpolation over triangular irregular networks (TIN);
polynomial;
splines;
kernel, e.g. Lanczos interpolation;
Thiessen or nearest neighbor;
inverse distance weighting (IDW), a generalization of nearest
neighbor interpolation;
kriging and related geostatistics methods;
geographically weighted regression (GWR);
...
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
138 / 164
Comparison of interpolation techniques
The figure below illustrates a set of distinct interpolation techniques in
one dimension. The dots represent discrete measurements and the
graphs represents the interpolated values. Most techniques illustrated
below preserve the observation at the original locations but lead to
different interpolated values for locations in between.
Thiessen
Splines
IDW
Kriging
Polynomial
Kriging
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
139 / 164
Some interpolation techniques for spatial data sets
1
Nearest neighbor: this technique of interpolation associates to
each point (x, y ) the measurement zi of the closest point (xi , yi );
Spatially, this corresponds to attribute the value zi to all locations
with the Voronoi polygon (see Slide 141 ) of (xi , yi ); The major
drawback of this approach is that it leads to a discontinuous
spatial interpolation;
2
Linear interpolation over an irregular triangular network; The
assumptions of are simple like for the previous technique; the
interpolation is continuous but not smooth (see Slide 141);
3
Moving window; This large family of techniques uses a moving
window (or kernel) over the area to compute the convolution of the
kernel with the available measurements; The shape and the
distribution of weights of the kernel define the output: if the kernel
is uniform the technique is called moving average; if the kernel
depends on the inverse of the distance to its center the technique
is called inverse distance weighting (see Slide 142).
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
140 / 164
Linear interpolation and triangular irregular networks (TIN)
Given a set of (xi , yi ) points (black dots in figure A), the Delaunay
triangulation fills their convex hull with triangles. Each triangle is
defined by three points such that their circumcircle doesn’t contain a
fourth point. The triangles’ circumcenters (red dots) define the
Voronoi polygons of the set of points (red polygons in B).
For each triangle with vertices with u, v , w indices, the linear
interpolation is given by the linear function z = a + b x + c y such that
zi = a + b xi + c yi for all i = u, v , w.
(A)
Manuel Campagnolo (ISA/ULisboa)
(B)
SIGs com R
5 a 10 de Fevereiro de 2015
141 / 164
Spatial interpolation with a moving window
A kernel k is a 2-dimensional distribution like any of the illustrated
below. Its support is the domain where k > 0. The z estimate at
location (x0 , y0 ) is givenR by
R the following integral which is called a
convolution of k and f :
k (x − x0 , y − y0 ) × f (x, y ) dxdy .
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
142 / 164
Spatial interpolation with a moving window (cont.)
Given a set of measurements {z1 , . . . , zn } at spatial locations
{(x1 , y1 ), . . . , (xn , yn )} consider the function f (x, y ) = zi when (x, y ) is
one of the locations and f = 0 otherwise.
For a given kernel k ≥ 0, at any location (x0 , y0 ) the convolution of k
with f is a weighted average of all measurements that fall inside the
support of k . The weights are given by the shape of k .
The kernel is “moved” across
the space to compute the
convolution at each location
(x, y ). The result is a spatial
interpolation of the available
measurements.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
143 / 164
Spatial Interpolation: an example of application
Ecology data sources for GIS are frequently derived by some kind of
spatial interpolation technique. Among available data sources, many
being freely-available, for inclusion in GIS, it is not always clear which
data are native (as for most of remote sensing data) and which data
are the result of spatial interpolation of point measurements.
Four types of data sources that provide GIS readable data.
Terrestrial – e.g. Atlas of the Biosphere, South And Central
America Environmental Data, MODIS-Derived Satellite Data;
Climatic – e.g. IPCC Data Distribution Centre, Climatic Research
Unit of the University Of East Anglia, WorldClim Global Climate
Data;
Elevation, altitude, bathymetric and topographic – e.g. SRTM 90m
Digital Elevation Data and other sources at Earth Explorer;
Marine – e.g. Permanent Service For Mean Sea Level, World
Oceans Atlas.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
144 / 164
Spatial Interpolation: an example of application
Ecology data sources for GIS are frequently derived by some kind of
spatial interpolation technique. Among available data sources, many
being freely-available, for inclusion in GIS, it is not always clear which
data are native (as for most of remote sensing data) and which data
are the result of spatial interpolation of point measurements.
Four types of data sources that provide GIS readable data.
Terrestrial – e.g. Atlas of the Biosphere, South And Central
America Environmental Data, MODIS-Derived Satellite Data;
Climatic – e.g. IPCC Data Distribution Centre, Climatic Research
Unit of the University Of East Anglia, WorldClim Global Climate
Data;
Elevation, altitude, bathymetric and topographic – e.g. SRTM 90m
Digital Elevation Data and other sources at Earth Explorer;
Marine – e.g. Permanent Service For Mean Sea Level, World
Oceans Atlas.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
144 / 164
Spatial Interpolation: the case of WorldClim
WorldClim is a climatic data set based on punctual temperature and
precipitation time series.
It is a raster data set with ≈1km nominal resolution (30 arc sec),
derived from a set of weather station data, which is widely used in
ecological modeling. The point data come from the following sources.
Global Historical Climate Network Dataset (GHCN): over than
20000 weather stations for precipitation and over 7000 for
temperature;
World Meteorological Organization climatological normals: over
4000 and 3000 stations;
FAOCLIM: over 27000 and 20000 stations;
Other regional databases.
Overall, and excluding duplicate records, Worldclim is based on over
47000 locations for precipitation, 25000 for mean temperature, and
14000 for min and max temperatures.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
145 / 164
Spatial Interpolation: some limitations of WorldClim
1
Although the set of weather stations is dense in parts of the world,
the spatial distributions of stations is very irregular.
Figure: Spatial distribution of weather stations of the Global Historical
Climate Network Dataset.
2
Worldclim is a raster dataset that covers all land areas – for most
locations it only provides estimates of the climatic variables
values.
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
146 / 164
Interpolação espacial com técnicas de Kriging
Sugestão de referência: http://www.css.cornell.edu/
faculty/dgr2/teach/R/R_ck.pdf.
Para ilustrar este tópicos vão ser usados os seguintes conjuntos de
dados localizados na Península Ibérica e obtidos através do Global
Historical Climate Network Dataset.
1
2
1
2
3
4
5
cdg de pontos stations_iberia_wgs84; os atributos jan,
fev, são as temperaturas médias (em 0.01o C) medidas em
estações do Global Historical Climate Network Dataset;
cdg de polígonos Iberian_Peninsula.
temps<−readOGR ( dsn=getwd ( ) , l a y e r = " s t a t i o n s _ i b e r i a _wgs84 " ,
encoding= " ISO8859−1" )
tem ps@pro j 4s t ri n g
i b e r i a<−readOGR ( dsn=getwd ( ) , l a y e r = " I b e r i a n P e n i n s u l a " , encoding= "
ISO8859−1" )
iberia@proj4string
tempib<−spTransform ( temps , i b e r i a @ p r o j 4 s t r i n g )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
147 / 164
Interpolação espacial com técnicas de Kriging
Eliminar observações que não são “continentais” com uma selecção
por localização:
1
2
3
4
5
6
7
8
9
10
11
p l o t ( tempib )
p l o t ( i b e r i a , add=TRUE)
# vamos s e l e c c i o n a r apenas as " c o n t i n e n t a i s "
i t s <− g I n t e r s e c t s ( tempib , i b e r i a , b y i d =TRUE) # devolve m a t r i z 2 * 56
# c a l c u l a r v e c t o r l ó g i c o dos pontos que estão em 1 dos 2
p o l í g o n o s " c o n t i n e n t a i s " ( P o r t u g a l , Espanha )
c t l <−as . l o g i c a l ( a p p l y ( i t s , 2 , max ) )
# s e l e c c i o n a r de tempib os pontos " c o n t i n e n t a i s " :
tempib<−tempib [ c t l , ] # tem 45 pontos
xy<−c o o r d i n a t e s ( tempib )
plot ( iberia )
t e x t ( xy [ , 1 ] , xy [ , 2 ] , tempib@data$name , cex = . 6 )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
148 / 164
Interpolação espacial com técnicas de Kriging
Pré-processar dados: criar uma grelha de pontos sobre a Península
Ibérica onde vai ser feita a interpolação espacial.
1
2
3
4
5
6
7
x y t<−data . frame ( c o o r d i n a t e s ( tempib ) )
colnames ( x y t )<−c ( " x " , " y " )
# Create sampling g r i d
x . range <− as . i n t e g e r ( range ( x y t [ , " x " ] ) )
y . range <− as . i n t e g e r ( range ( x y t [ , " y " ] ) )
grd <− expand . g r i d ( x = seq ( x . range [ 1 ] , x . range [ 2 ] , l e n g t h . o u t
=20) , y = seq ( y . range [ 1 ] , y . range [ 2 ] , l e n g t h . o u t =20) )
grd . sp <− S p a t i a l P o i n t s ( grd , p r o j 4 s t r i n g = i b e r i a @ p r o j 4 s t r i n g )
Carregar package gstat que contém as funções para Kriging.
Também contém funções para a técnica IDW mencionada
anteriormente.
1
l i b r a r y ( g s t a t ) # k r i g e , idw
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
149 / 164
Interpolação espacial com técnicas de Kriging
Oraganizar dados para Kriging. Para além da temperatura de Janeiro,
vamos usar a informação sobre a elevação:
1
2
e l e v<−as . numeric ( as . c h a r a c t e r ( tempib@data$ g r e l e v ) )
j a n<−as . numeric ( as . c h a r a c t e r ( tempib@data$ j a n ) )
Criar SpatialPointDataFrame com toda a informação necessária
para fazer Kriging:
1
2
3
4
5
6
7
8
# c r i a r S p a t i a l P o i n t s com essa informação
co<−data . frame ( x= x y t [ , " x " ] , y= x y t [ , " y " ] , j a n = jan , e l e v = e l e v ) #
45 * 2
# e l i m i n a r NAs
co<−co [ ! i s . na ( co$ e l e v ) , ]
# o comando abaixo c o n v e r t e a data . frame em
SpatialPointDataFrame
c o o r d i n a t e s ( co ) <− ~ x+y
# a s s o c i a r um CRS a co
p r o j 4 s t r i n g ( co ) <− i b e r i a @ p r o j 4 s t r i n g
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
150 / 164
Interpolação espacial com técnicas de Kriging
Primeiro modelo: modela-se unicamente a temperatura usando as
variáveis explicativas x e y , sem ter em conta a autocorrelação
espacial entre as observações da temperatura. Esta técnica é
equivalente a fazer uma regressão linear múltipla da temperatura
sobre x e y . O comando seguinte faz a predição das temperaturas
nos locais definidos em grid.sp:
1
k<− k r i g e ( j a n ~x+y , l o c a t i o n s =co , newdata=grd . sp )
Para mostrar o resultado apenas sobre a Península Ibérica pode
seleccionar-se as estimativas obtidas através de uma selecção por
localização:
1
2
3
4
5
6
# s e l e c c i o n a r os pontos de p r e d i ç ã o no i n t e r i o r da p e n í n s u l a
i t s <− g I n t e r s e c t s ( k , i b e r i a , b y i d =TRUE) # devolve m a t r i z 2 * 56
# c a l c u l a r v e c t o r l ó g i c o dos pontos que estão em 1 dos 2
p o l í g o n o s " c o n t i n e n t a i s " ( P o r t u g a l , Espanha )
i t s p i <−as . l o g i c a l ( a p p l y ( i t s , 2 , max ) )
# s e l e c c i o n a r de tempib os pontos " c o n t i n e n t a i s " :
k<−k [ i t s p i , ]
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
151 / 164
1
2
3
4
Interpolação espacial com técnicas de Kriging
Construir figura:
xyk<−c o o r d i n a t e s ( k )
plot ( iberia )
text ( xyt [ , " x " ] , xyt [ , " y " ] ,
round ( co@data [ , " j a n " ] )
, c o l = " b l u e " , cex = . 5 )
t e x t ( xyk [ , 1 ] , xyk [ , 2 ] , round
( k@data [ , " var1 . pred " ] )
, cex = . 5 )
O resultado é claramente
mau...
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
152 / 164
Interpolação espacial com técnicas de Kriging
Segundo modelo: modela-se unicamente a temperatura usando as
variáveis explicativas x e y , tendo em consideração a autocorrelação
espacial entre as observações da temperatura. O primeiro aspecto a
considerar é então a construção de um modelo para essa
autocorrelação espacial.
Construção do variograma:
1
2
3
4
# variograma , com pressuposto de a n i s o t r o p i a , l i m i t a d o a 1 / 3 da
d i s t â n c i a maior e 15 " l a g s "
p l o t ( variogram ( j a n ~x+y , l o c a t i o n s =co ) )
# variograma , com pressuposto de a n i s o t r o p i a , l i m i t a d o a 1 / 3 da
d i s t â n c i a maior e com l a g de 20km
p l o t ( variogram ( j a n ~x+y , l o c a t i o n s =co , w i d t h =20000) )
O passo seguinte é ajustar um modelo ao variograma e estimar os
parâmetros dess modelo. Em gstat existem uma séri de modelos
disponíveis, cuja lista é dada pelo comando:
1
vgm ( ) # g s t a t
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
153 / 164
Interpolação espacial com técnicas de Kriging
O variograma é parametrizado por sill, range e nugget. Na figura
também se ilustram alguns dos modelos mais utilizados:
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
154 / 164
1
2
3
4
5
Interpolação espacial com técnicas de Kriging
Ajustar um modelo ao variograma
com a função fit.variogram:
# a j u s t a r parâmetros do
variograma depois de
e s c o l h e r o modelo
v<−variogram ( j a n ~x+y , l o c a t i o n s =
co )
p l o t ( v ) # para a j u d a r a
e s c o l h e r ponto i n i c i a l da
pesquisa com vgm
v f<− f i t . variogram ( v , vgm ( p s i l l
=100000 , model= " Exp " , range
=300000 , nugget =0) )
plot (v , vf )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
155 / 164
Interpolação espacial com técnicas de Kriging
O comando seguinte faz a predição das temperaturas nos locais
definidos em grid.sp usando o modelo para o variograma ajustado
atrás vf:
1
k t<− k r i g e ( j a n ~x+y , l o c a t i o n s =co , newdata=grd . sp , model= v f )
Eliminar as predições fora da Península e construir gráfico:
1
2
3
4
5
6
7
8
i t s <− g I n t e r s e c t s ( k t , i b e r i a , b y i d =TRUE)
i t s p i <−as . l o g i c a l ( a p p l y ( i t s , 2 , max ) )
k t<− k t [ i t s p i , ] # tem 59 pontos dos 100 i n i c i a i s
x y k t<−c o o r d i n a t e s ( k t )
# imagem
plot ( iberia )
t e x t ( x y t [ , " x " ] , x y t [ , " y " ] , round ( co@data [ , " j a n " ] ) , c o l = " b l u e " , cex
=.5)
t e x t ( x y k t [ , 1 ] , x y k t [ , 2 ] , round ( kt@data [ , " var1 . pred " ] ) , cex = . 5 )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
156 / 164
Interpolação espacial com técnicas de Kriging
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
157 / 164
1
2
3
Interpolação espacial com técnicas de Kriging
Uma forma de fazer a validação
dos resultados é por validação
cruzada. No exemplo abaixo
parte-se o conjunto de pontos em
10 grupos para estimar os erros do
modelo.
# cv é S p a t i a l P o i n t s D a t a F r a m e
cv<− k r i g e . cv ( j a n ~x+y , l o c a t i o n s =
co , model= v f , n f o l d =10)
bubble ( cv , " r e s i d u a l " )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
158 / 164
Interpolação espacial com técnicas de Kriging
Terceiro modelo: modela-se a temperatura usando as variáveis
explicativas x e y e também a co-variável elevação, tendo em
consideração a autocorrelação espacial dessas variáveis e a
correlação entre temperatura e elevação. O primeiro aspecto a
considerar é então a construção de um modelo para esses três tipos
de correlação.
Construção do variograma para a temperatura (como atrás):
1
2
3
p l o t ( variogram ( j a n ~x+y , l o c a t i o n s =co ) )
v f t <− f i t . variogram ( v , vgm ( p s i l l =100000 , model= " Exp " , range =400000 ,
nugget =0) )
plot (v , vft )
Aplicar o mesmo procedimento para o variograma da elevação:
1
2
3
4
# a j u s t a r variograma
p l o t ( variogram ( e l e v ~x+y , l o c a t i o n s =co ) )
v f e<− f i t . variogram ( v , vgm ( p s i l l =200000 , model= " Exp " , range =400000 ,
nugget =0) )
p l o t ( v , vfe )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
159 / 164
Interpolação espacial com técnicas de Kriging
No package gstat é necessário construir uma estrutura única para
suportar a variável resposta e a(s) co-variável(eis):
1
2
3
g <− g s t a t ( NULL , i d = " j a n " , form = j a n ~ x+y , data=co )
g <− g s t a t ( g , i d = " e l e v " , form = e l e v ~ x+y , data=co )
g <− g s t a t ( g , i d = " j a n " , model = v f t , f i l l . a l l =T )
Em seguida podem ser ajustados os três modelos de correlação. O
caso mais simples consiste em ajustar o que é designado por modelo
linear de co-regionalização:
1
2
3
4
v . c r o s s <− variogram ( g )
p l o t ( v . cross , p l =T )
g <− f i t . lmc ( v . cross , g )
p l o t ( variogram ( g ) , model=g$model )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
160 / 164
Interpolação espacial com técnicas de Kriging
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
161 / 164
Interpolação espacial com técnicas de Kriging
O comando seguinte faz a predição das temperaturas usando a
co-variável elevação nos locais definidos em grid.sp. O modelo para
o variograma é o modelo linear de co-regionalização ajustado atrás.
1
kc <− p r e d i c t . g s t a t ( g , newdata=grd . sp )
Eliminar as predições fora da Península e construir gráfico:
1
2
3
4
5
6
7
8
i t s <− g I n t e r s e c t s ( kc , i b e r i a , b y i d =TRUE)
i t s p i <−as . l o g i c a l ( a p p l y ( i t s , 2 , max ) )
kc<−kc [ i t s p i , ]
xykc<−c o o r d i n a t e s ( kc )
# imagem
plot ( iberia )
t e x t ( x y t [ , " x " ] , x y t [ , " y " ] , round ( co@data [ , " j a n " ] ) , c o l = " b l u e " , cex
=.5)
t e x t ( xykc [ , 1 ] , xykc [ , 2 ] , round ( kc@data [ , " j a n . pred " ] ) , cex = . 5 )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
162 / 164
1
2
Interpolação espacial com técnicas de Kriging
Novamente, a validação dos
resultados é feita por validação
cruzada, particionando os dados
em 10 grupos:
cvc <− g s t a t . cv ( g , n f o l d =10)
bubble ( cvc , " r e s i d u a l " )
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
163 / 164
1
2
3
4
5
6
7
Interpolação espacial com técnicas de Kriging
Finalmente comparam-se os
resíduos obtidos usando ou não a
a u x t<−cv@data
co-variável temperatura. O modelo
auxc<−cvc@data
com co-variável tem menores
p l o t ( a u x t $observed , a u x t $ var1 .
resíduos.
pred , x l a b = " t e m p e r a tu r a
observada " , y l a b = "
t em p e r a t u r a estimada " , pch
=16 , asp =1)
p o i n t s ( auxc$observed , auxc$ j a n .
pred , c o l = " red " , pch =16)
abline (0 ,1)
t e x t ( min ( a u x t $observed ) , max (
a u x t $ var1 . pred ) , " sem co−
r e g i o n a l i z a ç ã o " , pos =4)
t e x t ( min ( a u x t $observed ) ,
q u a n t i l e ( a u x t $ var1 . pred
, 0 . 8 ) , " com co−
r e g i o n a l i z a ç ã o " , c o l = " red " ,
pos =4)
Manuel Campagnolo (ISA/ULisboa)
SIGs com R
5 a 10 de Fevereiro de 2015
164 / 164