How measure the mean distance of geospatial points in a map?

General
,
1

I don't understand spatial.data.

What I have: data.frame enterprises with the columns: id, parent_subsidiary, city_cod.

What I need: the mean and the max distance from the parent's city to the subsidiary cities.

Ex:

    id         |     mean_dist     | max_dist
 2152          |         25km      |     50km    
 3333          |        110km      |    180km  
 6666          |          0km      |      0km

Than, I need to do a map.

What I did : I have been trying to follow this: r - Geographic / geospatial distance between 2 lists of lat/lon points (coordinates) - Stack Overflow

library("tidyverse")
library("sf")
# library("brazilmaps")   not working anymore
library("geobr")


parent <- enterprises %>% filter(parent_subsidiary==1) 
subsidiary <- enterprises %>% filter(parent_subsidiary==2) 

# Cities - polygons 
m_city_br <- read_municipality(code_muni="all", year=2019)

# or shp_city<- st_read("/BR_Municipios_2019.shp")

# data.frame with the column geom
map_parent  <- left_join(parent, m_city_br, by=c("city_cod"="code_muni"))
map_subsidiary <- left_join(subsidiary, m_city_br, by=c("city_cod"="code_muni"))



st_distance(map_parent$geom[1],map_subsidiary$geom[2]) %>% units::set_units(km)
# it took a long time and the result is different from google.maps
# is it ok?!


# To do by ID -- I also stucked here

distance_p_s <- data.frame(id=as.numeric(),subsidiar=as.numeric(),mean_dist=as.numeric(),max_dist=as.numeric())

id_v <- as.vector(parent$id)



for (i in 1:length(id_v)){
  
 
  test_p <- map_parent %>% filter(id==id_v[i])  
  test_s <- map_subsidiary %>% filter(id==id_v[i])
  total <- 0
  value <- 0
  max <- 0
  l <- 0
  
  l <- nrow(test_s)

      for (j in 1:l){

         value <- as.numeric(round(st_distance(test_p$geom[1],test_s$geom[j]) %>% units::set_units(km),2))
          
         total <- total + value
         ifelse(value>max,max<-value,NA)
      }
  

  mean_dist <- total/l
  done <- data.frame(id=id[i],subsidiary=l,mean_dist=round(mean_dist,2),max_dist=max)
  distance_p_s <- rbind(distance_p_s,done)
  
  rm(done)
  
}
}

And I'm already stuck. I don't know how to proceed with 'longitude' and 'latitude' .
Can I calculate the centroid of the cities and than calculate the distance?

Data: structure(list(id = c("1111", "1111", "1111", "1111", "232", "232", "232", "232", "3123", "3123", "4455", "4455", "686", "333", "333", "14112", "14112", "14112", "3633", "3633","77172","77172"), parent_subsidiary = c("1","2", "2", "2", "1", "2", "2", "2", "1", "2", "1", "2", "1", "2", "1", "1", "2", "2", "1", "2","1","2"), city_cod = c(4305801L,4202404L, 4314803L, 4314902L, 4318705L, 1303403L, 4304507L, 4314100L, 2408102L, 3144409L, 5208707L, 4205407L, 5210000L, 3203908L, 3518800L, 3118601L, 4217303L, 3118601L, 5003702L, 5205109L,4111407L,4110102L)), row.names = c(NA, 22L), class = "data.frame")``

Thanks a lot !
PS: this is Brazilian cities

2

This package uses simple featuresand thesf` package has everything needed for finding centroids, plotting distances between them and can be combined with data frames for thematic mapping.

1 Like
3

Thanks!

Do you know where can I find some easy examples?

4 
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
# from the vignette
nc <- st_read(system.file("shape/nc.shp", package="sf"))
#> Reading layer `nc' from data source `/usr/local/lib/R/site-library/sf/shape/nc.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> geographic CRS: NAD27
x = st_transform(nc, 32119)
st_distance(x[c(1,4,22),], x[c(1, 33,55,56),])
#> Units: [m]
#>           [,1]     [,2]      [,3]     [,4]
#> [1,]      0.00 312184.9 128341.85 475623.3
#> [2,] 440561.15 114939.7 590434.80      0.0
#> [3,]  18944.03 352719.1  78756.89 517527.8

# https://journal.r-project.org/archive/2018/RJ-2018-009/RJ-2018-009.pdf
st_area(st_transform(nc[1, ], 2264)) %>% units::set_units(km^2)
#> 1137.598 [km^2]
 
# mine
nc_centroids <- st_centroid(nc)
#> Warning in st_centroid.sf(nc): st_centroid assumes attributes are constant over
#> geometries of x
#> Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
#> of_largest_polygon): st_centroid does not give correct centroids for longitude/
#> latitude data
x = st_transform(nc, 32119)
#alleghany <- nc_centroids[2]
st_distance(nc_centroids$geometry[1],nc_centroids$geometry[2])
#> Units: [m]
#>          [,1]
#> [1,] 34093.21

Created on 2020-08-19 by the reprex package (v0.3.0)

1 Like
5

Thanks a lot!

Now I just need to calculate the means and the max!

6

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.