4

Estoy intentando calcular la superficie que describen varias coordenadas (abajo descritas). Para ello uso el siguiente código

 library(sf)

  tbl <- readr::read_table2(
  Code Latitude Longitude
  C1    40.051302   -3.593308
  C2    40.043193   -3.599064
  H1    40.04933    -3.477574
  H2    40.085899   -3.214026
  H3    40.245852   -3.253943
  M1    40.03184    -3.34522
  M1    40.03184    -3.34522
  M2    40.050578   -3.50304
  M3    40.0459 -3.125
  M4    40.032298   -3.182801
  Z1    40.050052   -3.476076
  Z2    40.04472    -3.131
  L1    40.03185    -3.3453
  L2    40.048758   -3.476462
  L3    40.050352   -3.125794
  L4    40.245459   -3.252617
  Q1    39.59552    -3.39183
  Q2    40.225155   -3.530613
  Q3    40.320294   -3.205903
  Q4    40.345348   -3.137145
  E1    40.059138   -3.500323
  E2    40.234966   -3.536191
  E3    40.089637   -3.196057
  E4    40.494167   -3.131139") 

  tbl[, c(3,2)] %>%
    as.matrix() %>%
    list() %>%
    st_polygon() %>%
    st_sfc() %>% 
    st_set_crs(4326) %>% 
    st_area()

pero obtengo el siguiente error:

 Error in MtrxSet(x, dim, type = "POLYGON", needClosed = TRUE) : 
   polygons not (all) closed

¿Cómo se puede solucionar el error o que alternativas hay?

Mis puntos dibujan este área, me gustaría saber que superficie tienen.

introducir la descripción de la imagen aquí Gracias de antemano

3
  • La primera pregunta es sobre los datos mismos. Si le estás diciendo a MtrxSet que needClosed=TRUE, estás seguro que los datos de entrada contienen polígono(s) cerrados? Lo otro que puede estar pasando es algo así en donde los puntos inicial y final son diferentes por líos con el punto flotante.
    – Alfabravo
    el 7 sep. 2020 a las 13:07
  • Creo que la pregunta escapa al algoritmo como tal. ¿Cómo se calcula el área de un polígono que no está cerrado?
    – Alfabravo
    el 7 sep. 2020 a las 14:34
  • @Alfabravo el área que dibujan mis coordenadas se muestra en la edición de esta pregunta. Por favor, revisela. Gracias el 7 sep. 2020 a las 16:15

1 respuesta 1

3

Primero que nada, por los datos que muestras, no se tratarían de los puntos de un polígono sino de puntos individuales, por lo que tu primer problema es como encontrar un área que "encierre" todos los puntos. Lo que se me ocurre, es usar chull() para seleccionar aquellos puntos más externos, y sobre estos sí, construir el área o polígono.

puntos_totales <- tbl[, c(3,2)] 
puntos_area <- puntos_totales[chull(puntos), ]
puntos_area

     Longitude Latitude
[1,] -3.125000 40.04590
[2,] -3.391830 39.59552
[3,] -3.599064 40.04319
[4,] -3.536191 40.23497
[5,] -3.131139 40.49417

Sin embargo, todavía es necesario "cerrar" el área, el truco más fácil es replicar el primer punto, y ahora sí:

puntos_area %>% 
  rbind(puntos_area[1,]) -> puntos_area

puntos_area %>% 
  as.matrix() %>%
  list() %>%
  st_polygon() -> area

area %>% 
  st_sfc() %>% 
  st_set_crs(4326) %>% 
  st_area()

2303709063 [m^2]

La representación gráfica, sería:

plot(area)
points(puntos_totales %>% 
         as.matrix() %>%
         st_multipoint(), col="red")

introducir la descripción de la imagen aquí

2
  • Genial Patricio, funciona genial. Solo una objeción para futuras visitas de este post, cambia la línea del código: puntos_area <- puntos[chull(puntos), ] por puntos_area <- puntos_totales [chull(puntos_totales ), ] En definitiva, funciona genial, muchas gracias el 9 sep. 2020 a las 15:48
  • Absolutamente Adrián, corregido, gracias. el 9 sep. 2020 a las 15:53

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.