É 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.8

O 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.

x_{p} =\frac{\sum_{i=1}^n(\frac{1}{d_{i}^{2}}\times x_{i})}{\sum_{i=1}^n(\frac{1}{d_{i}^{2}})}

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) {
  # converte graus pra radiano
  rad <- pi/180
  # raio medio da terra no equador em km
  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 

  # distancia em km
  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    NA

Muito 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