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

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: https://stackoverflow.com/a/31668415/9591001

``````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

# 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,map_subsidiary\$geom) %>% 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,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

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

1 Like

Thanks!

Do you know where can I find some easy examples?

``````library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
# from the vignette
#> 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
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
st_distance(nc_centroids\$geometry,nc_centroids\$geometry)
#> Units: [m]
#>          [,1]
#> [1,] 34093.21
``````

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

1 Like

Thanks a lot!

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

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