Mapa temático da região de estudo

Incluindo a localização das estações meteorológicas

Introdução

A caracterização da região de estudo e das estações meteorológicas (EM) empregadas em pesquisas ou aplicações de meteorologia (ou áreas afins) é imprescindível para o melhor entendimento e interpretação dos resultados de uma análise observacional. Neste tutorial demonstra-se como produzir um mapa temático da região de estudo incluindo a localização das EM (que proveram os dados utilizados na pesquisa ou trabalho) e um atributo associado as EM (e.g. a temperatura média do ar climatológica da EM).

O mapa temático permitirá a visualização de informações relativas a:

• localização das EM
• delimitação dos estados brasileiros
• altitude do terreno
• algum atributo associado as EM

Pré-requisitos

Pacotes e funções necessárias.

## para instalar um pacote use
# install.packages("raster")
## limpando espaço de trabalho
rm(list = ls()) 
library(devtools)
## função load_pcks
source("../R/load_pcks.R")
## pacotes requeridos
pcks <- c("dplyr", "ggplot2","viridis", "lubridate", "scales", "raster", "viridis", "ggrepel")
load_pcks(pcks)
source_url('https://gist.githubusercontent.com/jdtatsch/5ec17e6624c97a7bbb767b07b46174ab/raw/47b38e36cc582ecc426dbc395edeb540db53caba/gg_bubble.R')
options(stringsAsFactors = TRUE)

Dados

Delimitação dos estados brasileiros

A delimitação¹ dos estados brasileiros pode ser obtida a partir de um arquivo shapefile fornecido pelo IBGE com os limites dos estados na escala de 1:250000. Abaixo mostra-se o procedimento para baixar e importar esses dados no R.

limites2015 <- "ftp://geoftp.ibge.gov.br/cartas_e_mapas/bases_cartograficas_continuas/bc250/versao2015/Shapefile/Limites_v2015_20160809.zip"
# nome e caminho para o arquivo que será baixado, altere se necessário
(zip_file <- paste0("../data/", basename(limites2015)))

[1] "../data/Limites_v2015_20160809.zip"

# baixando arquivo compactado
download.file(limites2015, destfile = zip_file)

# diretório para descompactar
extract_dir <- gsub("\\.zip", "", zip_file)
extract_dir

[1] "../data/Limites_v2015_20160809"

# descompactando arquivo
unzip(zip_file, exdir = extract_dir)
# lista dos shapefiles contidos no arquivo compactado
shapefiles_list <- list.files(extract_dir,
                              pattern = "shp$", 
                              recursive = TRUE, 
                              full.names = TRUE)
shapefiles_list

 [1] "../data/Limites_v2015_20160809/Limite_v2015_2016-08-03/LIM_Municipio_A.shp"                   
 [2] "../data/Limites_v2015_20160809/Limite_v2015_2016-08-03/LIM_Outras_Unid_Protegidas_A.shp"      
 [3] "../data/Limites_v2015_20160809/Limite_v2015_2016-08-03/LIM_Pais_A.shp"                        
 [4] "../data/Limites_v2015_20160809/Limite_v2015_2016-08-03/LIM_Terra_Indigena_A.shp"              
 [5] "../data/Limites_v2015_20160809/Limite_v2015_2016-08-03/LIM_Terra_Indigena_P.shp"              
 [6] "../data/Limites_v2015_20160809/Limite_v2015_2016-08-03/LIM_Terra_Publica_A.shp"               
 [7] "../data/Limites_v2015_20160809/Limite_v2015_2016-08-03/LIM_Unidade_Conservacao_Nao_Snuc_A.shp"
 [8] "../data/Limites_v2015_20160809/Limite_v2015_2016-08-03/LIM_Unidade_Federacao_A.shp"           
 [9] "../data/Limites_v2015_20160809/Limite_v2015_2016-08-03/LIM_Unidade_Protecao_Integral_A.shp"   
[10] "../data/Limites_v2015_20160809/Limite_v2015_2016-08-03/LIM_Unidade_Protecao_Integral_P.shp"   
[11] "../data/Limites_v2015_20160809/Limite_v2015_2016-08-03/LIM_Unidade_Uso_Sustentavel_A.shp"     
[12] "../data/Limites_v2015_20160809/Limite_v2015_2016-08-03/LIM_Unidade_Uso_Sustentavel_P.shp"     

Importando arquivo shapefile com delimitação dos estados.

# nome do arquivo shapefile dos estados
br_states_file <- grep(x = shapefiles_list, 
                       pattern = "Unidade_Federacao", 
                       value = TRUE)
# importa shape
br_states <- shapefile(br_states_file)

Z-dimension discarded

# projeta para lonlat
br_states <- spTransform(br_states, CRSobj = CRS("+proj=longlat +ellps=WGS84"))
plot(br_states, axes = TRUE)

Para selecionar a delimitação de alguns estados podemos fazer o seguinte:

sul <- br_states[br_states@data$NOME %in% c("Rio Grande do Sul", "Santa Catarina", "Paraná"), ]

Para uso dessa informação na função que gera o gráfico temático da região de estudo precisamos converter os dados armazenados na forma de um SpatialGridDataFrame para um dataframe, o que pode ser feito com a função fortify() do pacote ggplot2.

# 'fortificando' os dados (conversão para dataframe)
class(sul)

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

sul_df <- fortify(sul)

Regions defined for each Polygons

head(sul_df)

class(sul_df)

[1] "data.frame"

Altitude do terreno

A altitude da região de interesse pode ser obtida através da função getData() do pacote raster que faz o download do modelo digital de elevação do terreno da base de dados do

dem <- getData(name = "alt", download=TRUE, path = "../data", country = "BRA")
dem

class       : RasterLayer 
dimensions  : 4716, 5388, 25409808  (nrow, ncol, ncell)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : -74.1, -29.2, -33.9, 5.4  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 
data source : /home/hidrometeorologista/UFSM/orientacoes/iniciacaoCientifica/caroline/work/inmet_Caroline/data/BRA_msk_alt.grd 
names       : BRA_msk_alt 
values      : -127, 2665  (min, max)

[1] 0.008333333 0.008333333

O objeto dem precisa ser convertido de raster para um data.frame, análogo ao que foi feito com função fortify() no objeto sul, em que as colunas contém as coordenadas espaciais do raster dem (lon e lat) e os valores de altitude armazenados nas células do dem. Este dataframe é um dado de entrada da função criada para produção do mapa temático com a função gg_bubble(), mostrada a seguir.

class(dem)

[1] "RasterLayer"
attr(,"package")
[1] "raster"

# extent(dem)
# minha_extent <- extent(-50, -40, -20, -10)
# plot(dem); plot(minha_extent, add = TRUE)
dem_sul <- crop(dem, sul)
# mascará para os estados
dem_sul <- mask(dem, sul)
#plot(dem_sul); plot(sul, add = TRUE)
dem_sul_df <- data.frame(lon = xFromCell(dem_sul, cell = 1:ncell(dem_sul)),
                         lat = yFromCell(dem_sul, cell = 1:ncell(dem_sul)),
                         alt = values(dem_sul))
#summary(dem_sul_df)
# removendo linhas de dados faltantes
dem_sul_df <- dem_sul_df[complete.cases(dem_sul_df), ]
head(dem_sul_df)

class(dem_sul_df)

[1] "data.frame"

Eventualmente a região de estudo pode cobrir mais de um país. Nesse caso precisaremos baixar os dados de altitude do terreno para outros países² e juntá-los em um único raster. Esse procedimento é facilmente realizado com a função mosaic() do pacote raster. Veja ?mosaic para mais detalhes.

Informações das EM

As informações das estações meteorológicas devem incluir as coordenadas lon, lat e algum atributo de interesse, como por exemplo, a disponibilidade de dados, a temperatura média do ar, etc.

Nesse exemplo usaremos um dataframe preparado previamente e salvo em arquivo rds³.

info <- readRDS(file = "../output/info_sumary_tair_sul.rds")
info

# estrutura dos dados
str(info)

Classes ‘tbl_df’, ‘tbl’ and 'data.frame':   82 obs. of  17 variables:
 $ site    : chr  "A801" "A802" "A803" "A805" ...
 $ tmax_med: num  24.8 22.1 24.8 24.7 24.5 ...
 $ tmin_med: num  16.1 15.3 14.8 15.1 18 ...
 $ dtr_med : num  8.66 6.77 9.94 9.6 6.49 ...
 $ sdate   : Date, format: "2000-09-22" "2001-11-16" ...
 $ edate   : Date, format: "2015-12-31" "2015-12-31" ...
 $ period  : num  15.3 14.1 14.1 14.1 12.9 12.9 9.6 9.3 9.1 8.9 ...
 $ max_tair: num  40 38 40 36 39 34 41 40 39 37 ...
 $ min_tair: num  0 0 -2 -2 4 0 0 -3 -4 -2 ...
 $ missing : num  4.7 13.1 4.4 17.5 10.6 27.3 3.5 1.7 13 3.3 ...
 $ long_gap: num  305 9177 1110 9097 6639 ...
 $ sdate_lg: POSIXct, format: "2001-05-23 06:00:00" "2004-05-07 09:00:00" ...
 $ name    : chr  "PORTO ALEGRE" "RIO GRANDE" "SANTA MARIA" "SANTO AUGUSTO" ...
 $ state   : chr  "RS" "RS" "RS" "RS" ...
 $ lon     : num  -51.2 -52.1 -53.7 -53.8 -48.6 ...
 $ lat     : num  -30.1 -32 -29.7 -27.9 -27.6 ...
 $ alt     : num  46.97 2.46 95 550 1.8 ...

# nome das variáveis do data.frame (tabela de dados)
names(info)

 [1] "site"     "tmax_med" "tmin_med" "dtr_med"  "sdate"    "edate"   
 [7] "period"   "max_tair" "min_tair" "missing"  "long_gap" "sdate_lg"
[13] "name"     "state"    "lon"      "lat"      "alt"     

Os dados armazenados no arquivo RDS poderiam ser obtidos a partir de um arquivo texto contendo as informações das EM. Os dados dessa tabela poderiam ser importados como mostrado abaixo.

info <- read.table(file = "mydata.txt", 
                   header = FALSE,  # dados tem cabeçalho TRUE/FALSE
                   sep = " ",       # separador das colunas (",", "\t", " ")
                   na.strings = "-999.9") # string repreentando dados faltantes
# se os dados não tiverem cabeçalho (header = FALSE)
names(info) <- c("site", "tmax_med", "tmin_med", ...)
# salvando como RDS para importação mais rápida
saveRDS(info, "info_sumary_tair_sul.rds")

Gráfico com período de dados

Para demonstrar o uso a da função gg_bubble() faremos dois mapas temáticos. O primeiro representará a variação da média da temperatura mínima do ar para cada EM. No plano de fundo será mostrado a altitude do terreno, sobreposto pelo contorno dos estados.

# mapa da % de dados faltantes em cada EM
tmin_plot <- gg_bubble(data = info      # tabela de info das EM
                       ,z = "tmin_med"         # coluna da tabela data
                       ,colors_z = viridis     # paleta para variação de z
                       ,limites = sul_df       # dataframe com contorno da região (gerado de algum shapefile)
                       ,raster_bg = dem_sul_df # raster de background (altitude, decividade, etc)
                       ,colors_bg = gray.colors # paleta de cores do raster 
                       
                       ,z_legend = "Tmin (°C)" # texto para legenda
                       ,text_color = "red")     # cor do texto para o identificador ("site") das estações 
tmin_plot

No segundo mapa temático, alteramos o atributo das EM para o período de dados disponível de cada EM e não especificaremos os argumentos raster_bg e colors_bg. Pelo mapa resultante observa-se que os polígonos dos estados foram coloridos com uma mesma cor (definida pelo argumento color_fill com valor default burlywood3). O tamanho do label usado para indicar o código das EM foi aumentado para ilustrar o posicionamento otimizado dos labels.

# para período de dados
period_plot <-  gg_bubble(data = info
                          ,z = "period"
                          ,colors_z = viridis
                          ,limites = sul_df 
                          #,raster_bg = mdet_rs_df 
                          #,colors_bg = gray.colors
                          ,z_legend = "Período (anos)"
                          ,text_color = "black"
                          ,text_size = 4) 
period_plot

Pronto! Agora você tem um mapa temático de alta qualidade para caracterização de sua região de estudo.

---

1. A delimitação das regiões administrativas de qualquer país pode ser obtida através da função getData() do pacote raster que faz download da base de dados limites administrativos globais (GADM). Para obter os polígonos dos estados devemos definir o argumento level = 1. Esse argumento indica o nível de sub-divisão administrativa (level = 0 para o polígono do país, 1 para estados e 2 para municípios).↩

2. Repetimos a operação realizada com a função getData() alterando o argumento country conforme os países de interesse.↩

3. arquivo de armazenamento de dados em formato binário utilizado pelo R.↩