1

I would be thankful if you could help me with the following problem. I have a map for US counties and each county has its own characteristic, say, soil conductivity. Then I have coordinates of the source and signal receiver. What I need to do is to measure the distance from my source to the receiver and save the information about the distance and my signal went through for each county. For example, my source is in county A and my receiver is in county B, so I need to save the distance it moves through county A and the distance it moves through county B. Number of counties it goes through might be more than 2. I have built the map and put coordinates of my source and my receiver on it.

library(sf)
library(tmap)
county_map <- st_read("cty_1930_gis.shp")
qtm(county_map, fill = "cndctmn")
tm_shape(county_map) + tm_polygons(alpha = 0.2, legend.show = FALSE) + tm_shape(stations1) + tm_bubbles(col = "red", size=0.02) + tm_basemap(server = c('OpenStreetMap'))

conductivity for counties map

enter image description here

my source and my receiver on the map

enter image description here

What function should I use to calculate the distance my signal moves in each of counties, based on their borders I have on the map?

Ronak Shah
  • 377,200
  • 20
  • 156
  • 213
Daniel Alexsandrovich
  • 25
  • 3

2 Answers2

0

This is an sf question, not a tmap question.

This is the approach I recommend:

  1. Create spatial lines with sf from origin (source) to destination (receiver). Make sure to use a map projection (crs) that approximate distances well (so do not use mercator/wgs84).
  2. Use sf::st_intersection with two arguments: the lines you just created and the spatial polygons.
  3. The result already tells you through which counties the signal went.
  4. Use st::st_distance to measure the distance for each intersected piece.

Not sure if you are already familiar with sf. If not there are plenty of SO topics that cover the details of one of these steps.

Martijn Tennekes
  • 1,951
  • 13
  • 19
  • Thank you for the recommendation. Actually, I have already found a way to do it in Python, but it works really long, like 1 hour for 10k pairs of coordinates without parallelizing. Not sure if the R package can do it faster – Daniel Alexsandrovich Apr 07 '21 at 20:14
0

Here is the code that answers my question in Python in case someone would like to use it

import geopandas
from geopandas import GeoSeries
from shapely.geometry import LineString

data = geopandas.read_file('cty_1930_gis.shp')
data = data.to_crs(epsg=4326)
line = LineString([(lon_x, lat_x), (lon_y, lat_y)])
intersections = map_data.intersection(line)
intersections = intersections[~intersections.is_empty]
intersections = intersections.to_crs(epsg=3815)
lengths = intersections.length

and than write a loop of vector function over all points of ones' interest. (For sure, it's important to choose the correct epsg code)

Daniel Alexsandrovich
  • 25
  • 3