Interpolação pelo inverso do quadrado da distância
É comum quando temos um determinado valor distribuído espacialmente e queremos estimá-lo para um ponto específico. Existem inúmeras formas de se chegar nesta estimativa, mas quero mostrar apenas uma neste post. O objetivo é estimar o quanto choveu em Itapetininga-SP, a partir de dados de chuva de outras 6 cidades próximas. Utilizaremos para isso os dados das estações automáticas do INMET.
Primeiro, vamos importar e visualizar os dados que temos disponível.
if (!require("pacman")) install.packages("pacman")
pacman::p_load(readr, dplyr, leaflet)# importa o arquivo os dados de chuva
dados <- read_csv2(
"https://raw.githubusercontent.com/italocegatta/italocegatta.github.io_source/master/content/dados/chuva_inmet.csv"
)
dados# # A tibble: 6 x 4
# cidade lon lat p
# <chr> <dbl> <dbl> <dbl>
# 1 Sorocaba -47.6 -23.4 27
# 2 Itapeva -48.9 -24.0 33.4
# 3 Sao Miguel Arcanjo -48.2 -23.9 34.6
# 4 Avare -48.9 -23.1 18.2
# 5 Piracicaba -47.6 -22.7 30.8
# 6 Barra Bonita -48.6 -22.5 42.8O mapa a seguir mostra o total de chuva registrado pela estação meteorológica de cada cidade no dia 26/04/2017. Nosso objetivo é estimar o quanto choveu em Itapetininga utilizando a interpolação pelo inverso do quadrado da distância ou IDW (Inverse Distance Weighting).
leaflet(dados) %>%
addTiles() %>%
addMarkers(-48.0530600, -23.5916700) %>%
addCircleMarkers(
~lon, ~lat,
radius = ~p * 0.8,
label = ~as.character(p),
popup = ~cidade,
fillOpacity = 0.6,
labelOptions = labelOptions(
style = list("color" = "white"),
offset = c(5, -10),
noHide = TRUE,
textOnly = TRUE,
direction = "bottom"
)
)
A expressão que define o método é dada abaixo. Basicamente considera-se o valor de cada vizinho ponderado pelo inverso da distância entre ele e o ponto de interesse. Assim, vizinhos distantes contribuem com menos peso para o valor final que vizinhos mais próximos.
onde: xp = valor interpolado; xi = valor da i-ésimo ponto vizinho; di = distância entre o i-ésimo ponto de vizinho e o ponto de interesse.
Agora que já definimos o método, vamos começar os cálculos. O primeiro valor calculado será a distância entre os pontos. Utilizaremos a formula de Haversine que retorna a distâncias entre dois pontos de uma esfera a partir de suas latitudes e longitudes.
haversine <- function(lon1, lat1, lon2, lat2) {
rad <- pi/180
R <- 6378.1
dlon <- (lon2 - lon1) * rad
dlat <- (lat2 - lat1) * rad
a <- (sin(dlat/2))^2 +
cos(lat1 * rad) *
cos(lat2 * rad) *
(sin(dlon/2))^2
c <- 2 * atan2(sqrt(a), sqrt(1 - a))
d <- R * c
return(d)
}dist <- dados %>%
mutate(d_itape = haversine(lon, lat, -48.0530600, -23.5916700))
dist# # A tibble: 6 x 5
# cidade lon lat p d_itape
# <chr> <dbl> <dbl> <dbl> <dbl>
# 1 Sorocaba -47.6 -23.4 27 51.2
# 2 Itapeva -48.9 -24.0 33.4 95.3
# 3 Sao Miguel Arcanjo -48.2 -23.9 34.6 31.1
# 4 Avare -48.9 -23.1 18.2 106.
# 5 Piracicaba -47.6 -22.7 30.8 108.
# 6 Barra Bonita -48.6 -22.5 42.8 135.O cálculo do IDW é relativamente simples, basta reproduzir a expressão do método.
idw <- function(x, dist, na.rm = TRUE) {
s1 <- sum(x / dist^2, na.rm = na.rm)
s2 <- sum(1 / dist^2, na.rm = na.rm)
return(s1 / s2)
}dados_itape <- dist %>%
add_row(
.,
cidade = "Itapetininga",
lon = -48.0530600,
lat = -23.5916700,
p = round(idw(.$p, .$d_itape), 1)
)
dados_itape# # A tibble: 7 x 5
# cidade lon lat p d_itape
# <chr> <dbl> <dbl> <dbl> <dbl>
# 1 Sorocaba -47.6 -23.4 27 51.2
# 2 Itapeva -48.9 -24.0 33.4 95.3
# 3 Sao Miguel Arcanjo -48.2 -23.9 34.6 31.1
# 4 Avare -48.9 -23.1 18.2 106.
# 5 Piracicaba -47.6 -22.7 30.8 108.
# 6 Barra Bonita -48.6 -22.5 42.8 135.
# 7 Itapetininga -48.1 -23.6 32.1 NAMuito bom, agora vamos retornar ao mapa e adicionar o quanto choveu em Itapetininga de acordo com a interpolação por IDW.
Caso tenha alguma dúvida ou sugestão sobre o post, fique à vontade para fazer um comentário ou me contatar por E-mail.
sessioninfo::session_info(c("readr", "dplyr", "leaflet"))# - Session info ----------------------------------------------------------
# setting value
# version R version 3.5.3 (2019-03-11)
# os Windows 10 x64
# system x86_64, mingw32
# ui RTerm
# language (EN)
# collate Portuguese_Brazil.1252
# ctype Portuguese_Brazil.1252
# tz America/Sao_Paulo
# date 2019-08-25
#
# - Packages --------------------------------------------------------------
# package * version date lib source
# assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.5.3)
# backports 1.1.4 2019-04-10 [1] CRAN (R 3.5.3)
# base64enc 0.1-3 2015-07-28 [1] CRAN (R 3.5.0)
# BH 1.69.0-1 2019-01-07 [1] CRAN (R 3.5.2)
# cli 1.1.0 2019-03-19 [1] CRAN (R 3.5.3)
# clipr 0.7.0 2019-07-23 [1] CRAN (R 3.5.3)
# colorspace 1.4-1 2019-03-18 [1] CRAN (R 3.5.3)
# crayon 1.3.4 2017-09-16 [1] CRAN (R 3.5.1)
# crosstalk 1.0.0 2016-12-21 [1] CRAN (R 3.5.1)
# digest 0.6.20 2019-07-04 [1] CRAN (R 3.5.3)
# dplyr * 0.8.3 2019-07-04 [1] CRAN (R 3.5.3)
# ellipsis 0.2.0.1 2019-07-02 [1] CRAN (R 3.5.3)
# fansi 0.4.0 2018-10-05 [1] CRAN (R 3.5.1)
# ggplot2 3.2.1 2019-08-10 [1] CRAN (R 3.5.3)
# glue 1.3.1 2019-03-12 [1] CRAN (R 3.5.3)
# gridExtra 2.3 2017-09-09 [1] CRAN (R 3.5.1)
# gtable 0.3.0 2019-03-25 [1] CRAN (R 3.5.3)
# hms 0.5.0 2019-07-09 [1] CRAN (R 3.5.3)
# htmltools 0.3.6 2017-04-28 [1] CRAN (R 3.5.1)
# htmlwidgets 1.3 2018-09-30 [1] CRAN (R 3.5.1)
# httpuv 1.5.1 2019-04-05 [1] CRAN (R 3.5.3)
# jsonlite 1.6 2018-12-07 [1] CRAN (R 3.5.2)
# labeling 0.3 2014-08-23 [1] CRAN (R 3.5.0)
# later 0.8.0 2019-02-11 [1] CRAN (R 3.5.2)
# lattice 0.20-38 2018-11-04 [2] CRAN (R 3.5.3)
# lazyeval 0.2.2 2019-03-15 [1] CRAN (R 3.5.3)
# leaflet * 2.0.2 2018-08-27 [1] CRAN (R 3.5.1)
# magrittr 1.5 2014-11-22 [1] CRAN (R 3.5.1)
# markdown 1.1 2019-08-07 [1] CRAN (R 3.5.3)
# MASS 7.3-51.1 2018-11-01 [2] CRAN (R 3.5.3)
# Matrix 1.2-17 2019-03-22 [1] CRAN (R 3.5.3)
# mgcv 1.8-28 2019-03-21 [1] CRAN (R 3.5.3)
# mime 0.7 2019-06-11 [1] CRAN (R 3.5.3)
# munsell 0.5.0 2018-06-12 [1] CRAN (R 3.5.1)
# nlme 3.1-137 2018-04-07 [2] CRAN (R 3.5.3)
# pillar 1.4.2 2019-06-29 [1] CRAN (R 3.5.3)
# pkgconfig 2.0.2 2018-08-16 [1] CRAN (R 3.5.1)
# plogr 0.2.0 2018-03-25 [1] CRAN (R 3.5.1)
# plyr 1.8.4 2016-06-08 [1] CRAN (R 3.5.1)
# png 0.1-7 2013-12-03 [1] CRAN (R 3.5.0)
# promises 1.0.1 2018-04-13 [1] CRAN (R 3.5.1)
# purrr 0.3.2 2019-03-15 [1] CRAN (R 3.5.3)
# R6 2.4.0 2019-02-14 [1] CRAN (R 3.5.2)
# raster 2.9-23 2019-07-11 [1] CRAN (R 3.5.3)
# RColorBrewer 1.1-2 2014-12-07 [1] CRAN (R 3.5.0)
# Rcpp 1.0.2 2019-07-25 [1] CRAN (R 3.5.3)
# readr * 1.3.1 2018-12-21 [1] CRAN (R 3.5.2)
# reshape2 1.4.3 2017-12-11 [1] CRAN (R 3.5.1)
# rlang 0.4.0 2019-06-25 [1] CRAN (R 3.5.3)
# scales 1.0.0 2018-08-09 [1] CRAN (R 3.5.1)
# shiny 1.3.2 2019-04-22 [1] CRAN (R 3.5.3)
# sourcetools 0.1.7 2018-04-25 [1] CRAN (R 3.5.1)
# sp 1.3-1 2018-06-05 [1] CRAN (R 3.5.1)
# stringi 1.4.3 2019-03-12 [1] CRAN (R 3.5.3)
# stringr 1.4.0 2019-02-10 [1] CRAN (R 3.5.2)
# tibble 2.1.3 2019-06-06 [1] CRAN (R 3.5.3)
# tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.5.1)
# utf8 1.1.4 2018-05-24 [1] CRAN (R 3.5.1)
# vctrs 0.2.0 2019-07-05 [1] CRAN (R 3.5.3)
# viridis 0.5.1 2018-03-29 [1] CRAN (R 3.5.1)
# viridisLite 0.3.0 2018-02-01 [1] CRAN (R 3.5.1)
# withr 2.1.2 2018-03-15 [1] CRAN (R 3.5.1)
# xfun 0.8 2019-06-25 [1] CRAN (R 3.5.3)
# xtable 1.8-4 2019-04-21 [1] CRAN (R 3.5.3)
# yaml 2.2.0 2018-07-25 [1] CRAN (R 3.5.1)
# zeallot 0.1.0 2018-01-28 [1] CRAN (R 3.5.2)
#
# [1] C:/Users/Italo/Documents/R/win-library/3.5
# [2] C:/Program Files/R/R-3.5.3/library