Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • New version of -geo2xy- on SSC

    Thanks to Kit Baum, a new version of geo2xy is now available on SSC. If you already use geo2xy, you can update to the new version using
    Code:
    adoupdate geo2xy
    For a new installation, type in Stata's Command window:
    Code:
    ssc install geo2xy

    Examples in help files -- a new paradigm:

    The update includes an overhauled help system where each projection gets its own help file. These help files implement a new paradigm for presenting clickable examples, something that is particularly important for examples that generate graphs in Stata since multiple graph options make up for long command lines that need to be wrapped around for legibility. The model used by Stata in its official help files is suboptimal because the examples do not include line continuation indicators nor do they include a change of the delimiter that marks the end of a command. This means that if the code is copied, it cannot run without some surgery. With geo2xy, each example is completely self-contained and is run via a simple (click to run) link. The example is presented as do-file contents with a format that is similar to code blocks on this forum. Users can copy the example to a do-file and it is readily usable.

    This new model should be of interest to anyone who writes help files for Stata commands. Each (click to run) link runs a companion command that extracts the code as is from the help file itself, copies it to a temporary do-file, and then runs the temp do-file. The process is similar to selecting lines in Stata's do-file editor and running them by clicking the Do button. This means that you can edit your help file, develop your examples from within, reload the help file to visualize the changes and click to run the examples.

    The recent updates to rangestat and rangejoin (both from SSC) and the forthcoming rangerun command also implement the same model.


    Added projection: Lambert Conformal Conic Projection:

    geo2xy now includes a Lambert projection, contributed by Michael Stepner, which will be of particular interest to my fellow Canadians. It is Statistics Canada's recommended projection to use when generating maps of Canada's boundaries.

    While geo2xy help examples will work with any version of Stata 11 or higher, Stata 15 is now the lingua franca of Statalist so I'll use the new Stata 15 Sp command suite to illustrate this new map projection. So that everyone can play along, I'll use a shapefile that can be found on the following Natural Resources Canada site and downloaded via
    Code:
    http://ftp.geogratis.gc.ca/pub/nrcan_rncan/vector/atlas/base/7.5m_g_shp.zip
    You first need to unzip the archive (it includes a bunch of files that you do not need, the ones to look for start with "pvp"). Then you run the new spshape2dta command to create a Stata-format shapefile, which consists of two datasets, one for the attributes of each geographic unit and the other for the coordinates of the points that form the shape of these geographic units. The Sp suite includes the grmap command, a variant of spmap (Pisati 2004, from SSC). To create a map with this shapefile of Canada, you would:

    Code:
    * unzip the content and save to the current directory
    unzipfile 7.5m_g_shp.zip
    
    * convert to a Stata-format shapefile
    spshape2dta pvp, replace
    
    * load the attribute file and note the spset information
    use pvp
    spset
    
    * make a map!
    grmap
    The resulting map looks like
    Click image for larger version

Name:	canada.png
Views:	1
Size:	475.7 KB
ID:	1396973



    Now before you file a complaint with the Canadian government or go back to the internet to look for a better map of Canada, let me assure you that this is a perfectly acceptable shapefile and it is in fact the type of shapefile that you want to acquire because the coordinates are latitudes and longitudes. Before you can use it for analysis, you will have to tell Stata that these are not planar coordinates using
    Code:
    spset, modify coordsys(latlong)
    Note that this does not change the coordinates, it's just a flag that Stata uses to decide which distance formula to use: Euclidean if the coordinates are planar or great-circle distance (Haversine formula) if the coordinates are lat/lon. And if you redraw the map after modifying the spset coordinates to latlong, the map will look exactly the same. Trust me, this is a good thing, because you want to use great-circle distances with your spatial analysis. If you had a globe in front of you, latitudes are parallel to the equator. Longitudes are lines that go from pole to pole. A great-circle distance is the equivalent of pulling a string tight between two points on the globe. It represents the shortest path between two points along the surface of the globe.

    But what about making a usable map of Canada? That's where geo2xy comes to the rescue. Again using Stata 15 syntax, you can generate a map that uses the same projection that Google Maps uses with the following code:

    Code:
    * load the lat/lon coordinates and convert to planar
    use "pvp_shp.dta", clear
    geo2xy _Y _X, replace
    save "pvp_shp_xy.dta", replace
    
    * go back to the attributes dataset and modify the spset to point to the projected coordinates
    use "pvp.dta", clear
    spset, modify shpfile(pvp_shp_xy)
     
    * make a map!
    grmap
    and the map will look like

    Click image for larger version

Name:	canada_google.png
Views:	1
Size:	1.03 MB
ID:	1396974


    If you think this does not look right, go to Google Maps and zoom out. If you want to see what's going on, here's another take on the same map, this time using standard graph commands. The tissot option is used to generate circles, all circles have the exact same radius (as measured using great-circle distances) and are aligned on a longitude and latitude grid. Note that in order to create this map, all you need is the coordinates part of the shapefile.

    Code:
    use "pvp_shp.dta", clear
    geo2xy _Y _X , gen(ylat xlon) tissot
    
    // show the projection details and compute the plot's height
    return list
    local yheight = 6 * `r(aspect)'
    
    line ylat xlon if !mi(_ID), lwidth(vthin) lcolor(gray) cmissing(n) ///
    ||  ///
    line ylat xlon if mi(_ID), lwidth(vthin) lcolor(eltblue) cmissing(n) ///
        xsize(6) ysize(`yheight') ///
        ylabel(minmax, nogrid) yscale(off) ///
        xlabel(minmax, nogrid) xscale(off) ///
        plotregion(margin(small)) graphregion(margin(small)) ///
        legend(off) name(canada_mercator, replace)
    and the resulting map:

    Click image for larger version

Name:	canada_google_tissot.png
Views:	1
Size:	1.22 MB
ID:	1396975


    I hope that the above map shows the problem of doing spatial analysis using planar coordinates; the distances will be correctly calculated with respect to the map itself but they do not represent the correct geographic distances. Since maps are invariably distorted representations of the surface of the earth, it is best to obtain data located using geographic coordinates. When it comes time to make a map, you can easily apply a projection using geo2xy.

    Back to Canada, what are we supposed to do? Well... listen to the advice of Canadians in the know and use the projection they recommend. Again, with the tissot option to see the map distortions and thanks to Michael Stepner who contributed the code for the Lambert projection:

    Code:
    use "pvp_shp.dta", clear
    geo2xy _Y _X , gen(ylat xlon) projection(lambert_sphere, 49 77 0 -91.866667) tissot
    
    // show the projection details and compute the plot's height
    return list
    local yheight = 6 * `r(aspect)'
    
    line ylat xlon if !mi(_ID), lwidth(vthin) lcolor(gray) cmissing(n) ///
    ||  ///
    line ylat xlon if mi(_ID), lwidth(vthin) lcolor(eltblue) cmissing(n) ///
        xsize(6) ysize(`yheight') ///
        ylabel(minmax, nogrid) yscale(off) ///
        xlabel(minmax, nogrid) xscale(off) ///
        plotregion(margin(small)) graphregion(margin(small)) ///
        legend(off) name(canada_lambert, replace)
    and the better map of Canada:

    Click image for larger version

Name:	canada_lambert_tissot.png
Views:	1
Size:	822.9 KB
ID:	1396976



    Fun with Maps:

    geo2xy includes a "Fun with maps" help file. You can see how to replicate this Wikipedia map of the world, graticule and all, using a geo2xy's albers projection. There's also an example of how to create a composite map of the 48 conterminous US states with Alaska and Hawaii. I scooped myself on that one in this recent post.

    Enjoy!

  • #2
    Very helpful

    Comment

    Working...
    X