Add scale bar (in meters) to scatterplot of GPS coordinates in Stata

Question

I am graphing a series of GPS points using twoway scatter in Stata. This is just to visualize where the points are for some quick validation. I would like to add a scale bar for both the x and y axes to get a sense of how far apart the points actually are.

Here is some code to input a few coordinates and graph a scatterplot, specifying the axis ranges:

* Input example data
    clear
    input key lat lon
    1   -17.01639   36.11884
    2   -17.01669   36.11916
    3   -17.01692   36.11666
    4   -17.01575   36.11607
    5   -17.01551   36.11619
    6   -17.01543   36.11665
    7   -17.01581   36.11706
    8   -17.01665   36.11672
    9   -17.01543   36.11612
    10  -17.01724   36.11668
    end

* Install Geodist (we will need this later)
    cap which geodist
        if _rc == 111 ssc install geodist   // Installs geodist if necessary

* Create y/x axis ranges using min and max latitude and longitude
    sort lat
    local ymin = lat[1]
    local ymax = lat[_N]
    local ydelta = (`ymax' - `ymin')/3

    sort lon
    local xmin = lon[1]
    local xmax = lon[_N]
    local xdelta = (`xmax' - `xmin')/3

* Graph points
    twoway scatter lat lon,                                 ///
        ylabel(`ymin'(`ydelta')`ymax', angle(forty_five))   ///
        xlabel(`xmin'(`xdelta')`xmax')

You should get a lovely basic scatterplot. This is great, but I have no idea how far the points actually are from each other. Now, we can get a sense of the distance between the minimum and maximum latitude/longitude using the following:

* Store x and y ranges in locals
    geodist `ymin' `xmin' `ymax' `xmin'
    local yrange = `r(distance)'

    geodist `ymin' `xmin' `ymin' `xmax'
    local xrange = `r(distance)'

For this data, y spans 0.2km and x spans 0.329km.

Note: I know that on a spherical projection of the earth the distance between ymin and ymax may not map perfectly onto subsets of that distance (e.g. the distance from ymin to ymax/2 may not be precisely the same as half the distance between ymin and ymax). But for what I need this visualization for, and the fact that these data are usally clustered within 2km of each other, its accurate enough.

So, how can I create a scale bar for each of x and y to show what, for example, 100m looks like? In this example, the y scale would look different than the x scale, so maybe I need to find a way to change the y-ranges and x-ranges so that they can represent the same scale?

Answer 1

If you are going to create a map with geographic coordinates, you first need to convert these to planar (x,y) coordinates. With geo2xy (from SSC), you can use the same projection that Google Maps uses so the map you create will have exactly the same proportion that you would observe on Google Maps. Since you want coordinates that are based on distance units, I chose a true Mercator projection which creates (x,y) coordinates in meters. The projection is centered at mid-longitude by default which means that x-coordinates will be in meters from the center of the map. The y-coordinates are in meters from the equator. To center them on the middle of the map, all you need is to offset them from the mean.

When creating a map in Stata, you need to be sure to create a map region that has the correct proportions (xsize() and ysize() are proportional to the distance).

* Input example data
clear
input key lat lon
1   -17.01639   36.11884
2   -17.01669   36.11916
3   -17.01692   36.11666
4   -17.01575   36.11607
5   -17.01551   36.11619
6   -17.01543   36.11665
7   -17.01581   36.11706
8   -17.01665   36.11672
9   -17.01543   36.11612
10  -17.01724   36.11668
end

* project lat/lon to (x,y) in meters using a Mercator projection 
geo2xy lat lon, gen(ylat xlon) project(mercator)

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

* the projection is in meters, centered on mid-longitude, center y coor as well
sum ylat, meanonly
replace ylat = ylat - r(mean)

scatter ylat xlon,  ///
        xsize(6) ysize(`yheight') ///
        ylabel(minmax, nogrid) yscale(off) ///
        xlabel(minmax, nogrid) xscale(off) ///
        plotregion(margin(small)) graphregion(margin(small)) ///
        legend(off) name(mercator_us, replace)

graph export my_map.png, width(600) replace

Here's the resulting map:

If you want to display axis units (in meters), you can remove the xlabel() and ylabel() lines. If you do so, Stata will likely rescale each axis to fit nice intervals and will therefore change the proportions a bit. Here's what the map looks like if I let Stata decide on the labels: