spmap guidance: is it possible to draw circles around dot

[Luca Aguzzoni]

I am trying to use spmap to produce a map on which I can draw a point and a circle of 100 km distance around that point.

Can spmap accomodate such a request?
Would I manually need to define polygons around my point of interest or is there a built in feature that could help me?

In case I have to define polygons, is there any command on stata that could do so?

many thanks,

Luca

[Robert Picard]

Here's a program that will calculate the coordinates of a destination point given distance and bearing from a starting point (source for the formula: http://www.movable-type.co.uk/scripts/latlong.html)

program define destpoint
    
    args dlat1 dlon1 d b lat2 lon2
    
    local R 6371
    
    tempname lat1 lon1 brng
    gen double `lat1' = `dlat1' * _pi / 180
    gen double `lon1' = `dlon1' * _pi / 180
    gen double `brng' = `b' * _pi / 180
    
    gen double `lat2' = asin(sin(`lat1') * cos(`d'/`R') + ///
                        cos(`lat1') * sin(`d'/`R') * cos(`brng'))
                        
    gen double `lon2' = `lon1' + atan2(sin(`brng') * sin(`d'/`R') * cos(`lat1'), ///
                                        cos(`d'/`R') - sin(`lat1') * sin(`lat2'))
                                        
    qui replace `lat2' = `lat2' * 180 / _pi
    qui replace `lon2' = `lon2' * 180 / _pi
    
end

You need to save the program in a file called "destpoint.ado" in Stata's current directory (or in your PLUS directory, see help sysdir) so that Stata will find it.

The following code uses the destpoint program defined above to create 360 points that are 100km from a center point (in Washington, DC). I also use geodist (available from SSC) to confirm that each point is 100km from the center point. Finally, the destination points are converted to the coordinates data file of a shapefile.

* Create 360 neighbor points around a center point (clat clon) at a
* distance of 100km
    
    clear
    set obs 360
    gen nid = _n
    
    local dtarget 100
    
    gen double clat = 38.892535
    gen double clon = -77.020730
    gen double bearing = _n - 1
    
    destpoint clat clon `dtarget' bearing nlat nlon
    
    
* Double-check that each point is at 100km from the center point
    
    geodist clat clon nlat nlon, gen(d) sphere
    assert abs(`dtarget' - d) < 1e-10
    drop d


* Show the points on a map

    scatter nlat nlon,  msize(vtiny) || scatter clat clon, ///
        aspect(1) msize(tiny) scale(.5) xlabel(minmax) ylabel(minmax)


* Convert these points into a polygon, the first obs in a shapefile is missing;
* the last obs loops back to the first non-missing point

    expand 2 if _n == 1 | _n == _N
    sort nid
    replace nlat = nlat[1] in l
    replace nlon = nlon[1] in l
    replace nlat = . in 1
    replace nlon = . in 1
    
    keep nlat nlon
    rename nlat _Y
    rename nlon _X
    gen _ID = 1
    order _ID _X _X
    
    
* Plot the polygon

    line _Y _X, aspect(1) xlabel(minmax) ylabel(minmax)

[Luca Aguzzoni]

Many Thanks Robert, your code worked perfectly fine, cheers, Luca

[Luca Aguzzoni]

Many Thanks Robert, your code worked perfectly fine, cheers, Luca

[Robert Picard]

Thanks for the closure and see the announcement for geocircles, a new command that evolved from my proposed solution in #2 and the new geo2xy command that converts latitude and longitude to cartesian coordinates so that geographic circles can be mapped with minimal distorsions. Both are now available on SSC.

[George Ford]

After trying a few other options, I wrote a simple program that worked for me. You can see it here: https://papers.ssrn.com/sol3/papers....act_id=3156061

[Robert Picard]

I had a look at what worked for you and its not directly related to the issue of this thread. Instead, you are trying to generate random geographical locations (lat/lon) within a certain distance of a central location. The program you created to perform this task generates random locations within the desired distance but these locations are not uniformly distributed over the circle's area. If you make a histogram of the latitudes, you will see that each bin contains about the same number of points which is not surprising given that you determine the latitude using

runiform(latmin,latmax)

where the interval bounds represent the latitudes of the 2 points directly above and below at the maximum distance allowed.

This is an interesting problem and I must admit that it has been a bit challenging to find a "correct" solution on the internets. Here's my implementation of a solution proposed by http://www.geomidpoint.com/random/calculation.html

* method http://www.geomidpoint.com/random/calculation.html
clear all
program georcircular
    args lat1 lon1 maxdist
    
    local R = 6371 / 1.609344        // earth's radius in miles
    
    // convert to radians
    local lat1 =  `lat1' * _pi / 180
    local lon1 =  `lon1' * _pi / 180
    local maxdist = `maxdist' / `R'
    
    tempvar rdist
    gen `rdist' = acos(runiform() * (cos(`maxdist') - 1) + 1)
    
    tempvar bearing
    gen `bearing' = 2 * c(pi) * runiform()
    
    gen double lat = asin(sin(`lat1') * cos(`rdist') + ///
        cos(`lat1') * sin(`rdist') * cos(`bearing'))
        
    gen double lon = `lon1' + atan2( sin(`bearing') * sin(`rdist') * ///
        cos(`lat1'), cos(`rdist')-sin(`lat1') * sin(lat))
        
    replace lat = lat * 180 / _pi
    replace lon = lon * 180 / _pi
    replace lon = lon - 360 if lon > 180
    replace lon = lon + 360 if lon < -180
    
end

set obs 100000
georcircular 33.520080 -86.801175 5

geodist 33.520080 -86.801175 lat lon, gen(d) sphere miles
sum d

geo2xy lat lon, gen(ylat xlon)
return list
local yheight = 6 * `r(aspect)'
scatter ylat xlon, msymbol(point) ///
                xsize(6) ysize(`yheight') ///
                ylabel(minmax, nogrid) yscale(off) ///
                xlabel(minmax, nogrid) xscale(off) ///
                plotregion(margin(small)) graphregion(margin(small)) ///
                legend(off)

graph export rcircular.png, width(1000) replace

and here are the results from the distance calculation (geodist is from SSC)

. geodist 33.520080 -86.801175 lat lon, gen(d) sphere miles

. sum d

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
           d |    100,000    3.333456    1.176165   .0181738   4.999983

.

Before you can create a map of geographic coordinates, you need to convert lat/lon to planar coordinates. You can do this with geo2xy (from SSC). The default projection is the same used by Google Maps. To create a map in Stata, you have to control the horizontal and vertical dimension of the graph to make sure that it has the correct proportions. The following confirms what we would expect (a square map and thus an aspect ratio of 1):

. return list

macros:
             r(aspect) : "1.000147550531355"
              r(model) : "Unit Sphere"
              r(pname) : "Web Mercator"
               r(zoom) : "0"

and finally, here's what the map of these 100K locations looks like:

[George Ford]

Much better. Thanks. Here's the overlay between the two (blue, mine) and (orange, Robert's).

[Jay Hornberg]

Hello,

Thanks to you both for all the work that you lay out in this thread.

I have a list of lat/long points that I need to alter for anonymization purposes.
I would like to replace each point with another randomly-generated point within a radius .3 km around the origin.

I am trying to incorporate the work you explain above, but am running into some roadblocks. Could you offer some guidance?

Many thanks,
Jay

Data example:
input id alat alon
1 42.103 -80.065
2 42.103 -80.065001
3 39.739 -75.605
4 37.499 -77.470
5 39.464 -77.962
6 27.844 -82.796
7 39.138 -84.509
8 38.271 -85.694
9 36.143 -86.805
10 35.832 -86.073
end

[Robert Picard]

You can use the destpoint program from #2 to calculate coordinates using a random bearing and distance from the original point. I'm using a distance of 50km so that the points do not show one on top of each other on the final map.

* Example generated by -dataex-. To install: ssc install dataex
clear
input float(id alat alon)
 1 42.103 -80.065
 2 42.103 -80.065
 3 39.739 -75.605
 4 37.499  -77.47
 5 39.464 -77.962
 6 27.844 -82.796
 7 39.138 -84.509
 8 38.271 -85.694
 9 36.143 -86.805
10 35.832 -86.073
end

* compute random bearing and distance (use 0-50km to create an example map)
set seed 31234
gen bearing = runiform(0,360)
gen dtarget = runiform(0,50)

* compute the new coordinates
destpoint alat alon dtarget bearing nlat nlon

* double-check that the new coordinates are at the correct distance
geodist alat alon nlat nlon, gen(d) sphere
gen check = abs(dtarget - d)

list

* make a map, use a Google Maps projection (ssc install geo2xy)
geo2xy alat alon, gen(yalat xalon)
geo2xy nlat nlon, gen(ynlat xnlon)
return list
local yheight = 6 * `r(aspect)'
scatter yalat xalon, msymbol(+) ///
|| ///
scatter ynlat xnlon, msymbol(x) ///
    xsize(6) ysize(`yheight') ///
    ylabel(minmax, nogrid) yscale(off) ///
    xlabel(minmax, nogrid) xscale(off) ///
    plotregion(margin(small)) graphregion(margin(small)) ///
    legend(off)

graph export "anonymization.png", width(1000) replace

The results from the list command:

. list

     +---------------------------------------------------------------------------------------------+
     | id     alat      alon    bearing    dtarget        nlat         nlon           d      check |
     |---------------------------------------------------------------------------------------------|
  1. |  1   42.103   -80.065   165.2247   12.06249   41.998101   -80.027776   12.062487   1.34e-12 |
  2. |  2   42.103   -80.065   276.8026   41.90966   42.146539   -80.569771    41.90966   6.47e-13 |
  3. |  3   39.739   -75.605    284.247   31.72536   39.808659   -75.964992   31.725365   3.87e-13 |
  4. |  4   37.499    -77.47   314.9474   4.569766   37.528028   -77.506678   4.5697656   3.26e-13 |
  5. |  5   39.464   -77.962   289.6342   34.27682   39.566973   -78.338628   34.276817   1.19e-12 |
     |---------------------------------------------------------------------------------------------|
  6. |  6   27.844   -82.796   152.8325   32.24474   27.585928   -82.646611   32.244743   3.77e-13 |
  7. |  7   39.138   -84.509   265.2328   37.81234   39.108925   -84.945732    37.81234   6.11e-13 |
  8. |  8   38.271   -85.694   310.5642    44.1819   38.528758    -86.07985     44.1819   2.56e-13 |
  9. |  9   36.143   -86.805    296.119   48.98247   36.335929   -87.295998   48.982468   4.62e-13 |
 10. | 10   35.832   -86.073   56.50816   37.75438   36.018854   -85.722911   37.754379   3.84e-13 |
     +---------------------------------------------------------------------------------------------+

and the map:

[Marlene Fruehling]

Hi all,
Thank you for all these careful discussion!
I am new to smap and am trying to map villages (level=0) and cities (level=1) on a map and connect the points with a line in case the two levels share the same ID (meaning that a person goes both to the village and the city). In the example below, I would want a line connecting the 4th and 35th observation, for example. There could be a point that needs to be connected with several others.
The data looks like this:
+--------------------------------------+
| level latitude longitude ID |
|--------------------------------------|
1. | 0 14.8984 -16.5432 102 |
2. | 0 14.9084 -16.4556 103 |
3. | 0 14.7076 -16.0167 103 |
4. | 0 15.3681 -16.1361 102 |

34. | 1 14.941592 -16.490471 101 |
35. | 1 14.941592 -16.490471 102 |


So far, I am merging this data to data generated from a shapefile and plot the points using smap.
use corr.dta, clear
merge 1:m id using data", nogen

spmap Shape_Area using coord0.dta, id(id) fcolor(white) point(xcoord(longitude) ycoord(latitude) by(level) fcolor(Blues) legenda(on)) legend(size(medium) position(1))


Now, as I said, I am struggling with getting the lines to connect observations that share the same ID.

Do you have any advice on how to do this?Thank you very much in advance,
Best,
Marlene