Geonear or geodist in a vector space (using direction as an input)

[Prachi Singh]

Hi there!

I am working at a problem where I have to calculate whether a sampled cluster falls in the downwind direction from a pollution source, that is the pollution from the source flows towards the sampled cluster. For this I will have to use wind direction along with cluster and pollution source's location. By using geonear / geodist command I can ascertain the distance between the polluting source and sampled cluster but I wouldn't know whether a sampled cluster falls in the downwind direction. So I wanted to know whether there is command which can use direction as an input as well while calculating distance in stata.

Thanks,
Prachi

[Robert Picard]

I'm a bit reluctant to give advice because the vagueness on what should be considered downwind from the pollution source. Presumably, time affects the progression and expansion of pollution plumes so intuitively, a solution would most likely require forming polygons that define the borders of each plume. You would then use geoinpoly (from SSC) to identify points of interest that fall within these plume polygons. This sounds pretty complicated to achieve and I have no suggestion on how to go about creating such polygons.

The concept of direction is a also a bit difficult to grasp when it relates to geographic coordinates. You can follow a line of constant bearing (using a compass) but that's not the shortest path unless you are headed due North or South. Conversely, if you travel in a straight line using the shortest path between two points, your bearing will constantly change as you make progress towards your destination. You can sidestep the issue by converting geographic coordinates (latitude and longitude) to planar coordinates using geo2xy (from SSC). If you use the default projection (Web Mercator), you will be working with a map where lines of constant bearing are represented as straight lines. You can then use basic trig to calculate distances and slopes.

Before trying to formulate a solution using geographic coordinates, I will create two demonstration datasets. The first includes two "polluter" locations and a wind direction (in degrees, north is 0) for each. The second contains random points around these two "polluters".

set seed 312
* create pollution sources
* Example generated by -dataex-. To install: ssc install dataex
clear
input byte pid str16 pname double(plat plon) int wind
1 "The Cube"         42.275892 -83.741885  90
2 "Lurie Bell Tower" 42.291988 -83.71621  180
end
save "pollution_sources.dta", replace

* create points around pollution sources
collapse plat plon
expand 5000
gen lat = plat + runiform(-.04, .04)
gen lon = plon + runiform(-.05, .05)
keep lat lon
gen id = _n
save "points.dta", replace

Neither geodist nor geonear (both from SSC) calculate bearings between points but the formula for the initial bearing can be found on this excellent web site by Chris Veness. The first step to a solution is to form all pairwise combinations of polluters with points of interests and then calculate the distance and initial bearing. From there, you can decide that points that are within a certain range of bearings around the wind direction (and within a certain distance) fall within the plume. To illustrate what's going on, I project the coordinates and create a map. Note that I kept things simple and did not address how to define a range of bearings that cross over the 360/0 degrees.

clear all

* formula from http://www.movable-type.co.uk/scripts/latlong.html
program initial_bearing
    
    args lat1 lon1 lat2 lon2 newvar
    
    tempname d2r r2d
    scalar `d2r' = _pi / 180
    scalar `r2d' = 180 / _pi

    gen `newvar' = atan2(sin((`lon2'-`lon1') * `d2r') * cos(`lat2' * `d2r') , ///
                            cos(`lat1' * `d2r') * sin(`lat2' * `d2r') - ///
                            sin(`lat1' * `d2r') * cos(`lat2' * `d2r') * ///
                            cos((`lon2'-`lon1') * `d2r'))
    
    // normalize atan2 results (-pi to pi) to range from 0 to 360 degrees
    replace `newvar' = mod((`newvar' * `r2d') + 360,3 60)
end

* form all combinations of sources and points
use "pollution_sources.dta", clear
cross using "points.dta"

* calculate distances and initial bearing
geodist plat plon lat lon, gen(d)
initial_bearing plat plon lat lon ib

* identify points within a certain range from the wind direction
gen inplume = inrange(ib, wind-10, wind+10) & d < 4

* project coordinates before creating a map
geo2xy plat plon, gen(yplat xplon)
geo2xy lat lon, gen(ylat xlon)

* show the projection details and compute the graph's height
return list
local yheight = 6 * `r(aspect)'

scatter yplat xplon, msymbol(O) mcolor(red) ///
|| ///
scatter ylat xlon if !inplume, msymbol(x) mcolor(green) ///
|| ///
scatter ylat xlon if inplume, msymbol(x) mcolor(red) ///
xsize(6) ysize(`yheight') ///
    ylabel(minmax, nogrid) yscale(off) ///
    xlabel(minmax, nogrid) xscale(off) ///
    plotregion(margin(small)) graphregion(margin(small)) ///
    legend(off)

graph export "plume.png", width(800) replace

and here's the resulting map:

[Prachi Singh]

Thank you so much Robert, I will go through your code and get back to you.

[Prachi Singh]

Dear Robert, This is exactly what I wanted to implement. Also I understand your reservations about initial bearing being different from final bearing which is especially true if the distance between the two points is large (for example Baghdad and Osaka in Chris's example) but in my case the points are much closer (maximum 100kms apart) so it shouldn't be as big a problem. Once again thanks for such a detailed answer!

[Anannya Gupta]

I have typed the exact code given by Robert. For some reason idk what, I have an error that says *__000008 invalid name
Please help.

[Andrew Musau]

The code works as it should do for me. Are you running the code in parts? You should run the entire code at once in a do-file. If the error persists, tell us about your version of Stata and copy and paste the exact output from the Results window. You may also want to

set trace on

to see exactly where the code fails.

[Anannya Gupta]

Hey, the error still persists.
ps: (am two hours old on this forum, I am still reading and trying to grasp the rules and FAQs of posting code and results, so apologies for my unintended oversight)

Stata version: 15.0

input byte pid str16 pname double(plat plon) int wind

pid pname plat plon wind
1. 1 carvajal 4.600 -74.15 170.4286
2. end

. save "pollution_sources.dta", replace
file pollution_sources.dta saved

. collapse plat plon

. expand 5000
(4,999 observations created)

. gen lat = plat + runiform(-0.04, 0.04)

. gen lon = plon + runiform(-0.05, 0.05)

. keep lat lon

. gen id = _n

. save "points.dta", replace
file points.dta saved

. clear all

. args lat1 lon1 lat2 lon2 newvar

. tempname d2r r2d

. scalar `d2r' = _pi / 180

. scalar `r2d' = 180 / _pi

. gen `newvar' = atan2(sin((`lon2'-`lon1') * `d2r') * cos(`lat2' * `d2r') , cos(`lat1' * `d2r') * sin(`lat2' * `d2r') - sin(`lat1' * `d2r') * cos(`lat
> 2' * `d2r') * cos((`lon2'-`lon1') * `d2r'))
*__000000 invalid name
r(198);

end of do-file

r(198);

.

[Andrew Musau]

You have clearly edited the code in #2 to exclude

* formula from http://www.movable-type.co.uk/scripts/latlong.html
program initial_bearing

The comment is commented out and will not interfere with the running of the code, but the program command is critical. Run the code exactly as it is.