{smcl} {* *! version 1.0.0 31jan2017}{...} {cmd:help geo2xy_picard} {hline} {title:Title} {phang} {cmd:geo2xy} {hline 2} Convert latitude and longitude to cartesian (x,y) coordinates {title:Map projection} {phang} Picard's equidistant cylindrical projection {title:Syntax} {p 8 16 2} {cmd:geo2xy} {it:lat_var lon_var} {ifin} {cmd:,} {c -(}{opt gen:erate(y_lat x_lon)}{c |}{opt replace}{c )-} [ {opt proj:ection}{cmd:(}{opt picard} [,{help geo2xy##proj_name:proj_opts}]{cmd:)} {opt ti:ssot} ] {synoptset 42 tabbed}{...} {synopthdr} {synoptline} {p2coldent :* {opt gen:erate(y_lat x_lon)}}New variable names for the projected coordinates (note the order!){p_end} {p2coldent :* {opt replace}}Replace the values in {it:lat_var lon_var} with the projected coordinates{p_end} {synopt :{opt proj:ection(picard [,proj_opts])}}Specifies the projection and, optionally, projection parameters{p_end} {synopt :{opt ti:ssot}}Append Tissot's indicatrices (grid of circles){p_end} {synoptline} {pstd}* Either {opt gen:erate(y_lat x_lon)} or {opt replace} is required. {p2colreset} {synoptset 15 tabbed}{...} {synopthdr :proj_opts} {synoptline} {synopt :{it:lat1}}standard parallel (default is mid-latitude){p_end} {synoptline} {p2colreset} {title:Description} {pstd} On a unit sphere, a distance along a parallel is proportional to cos(latitude). So for several years, I have been making maps in Stata by plotting raw geographic latitudes and longitudes and adjusting the aspect ratio of the map so that the proportions are correct at mid-latitude. The map ratio is calculated as (maxlat-minlat) / (cos(midlat) * (maxlon-minlon)) {pstd} I did not quite understand that this qualifies as a projection until I started researching projections for {bf:geo2xy} and found out that this is exactly the same transformation that occurs with the {it:{help geo2xy_equidistant_cylindrical:Equidistant Cylindrical projection}}. {pstd} The simplicity of being able to plot raw geographic latitudes and longitudes while maintaining the correct map aspect ratio is appealing so I rolled it into {bf:geo2xy}. All you need to do is apply the aspect ratio returned by {bf:geo2xy}. See the example below. {title:Spheroid and (x,y) coordinates units} {pstd} This projection assumes that the geographic latitude and longitude describe locations on a sphere. The (x,y) coordinates are not projected and are simply clones of the original latitudes and longitudes. {title:Example} {pstd} This example requires {cmd:geo2xy}'s ancillary datasets in the current directory. Click {stata `"net get geo2xy, from("http://fmwww.bc.edu/repec/bocode/g")"':here} to download them. {pstd} In the following example, a crosshair is plotted using the coordinates of the Washington Monument in Washington, DC. By default, the standard parallel is set to mid-latitude but in this case, it is aligned with Washington's latitude. That makes the map's horizontal distances accurate at that latitude. As with the {it:{help geo2xy_equidistant_cylindrical:Equidistant Cylindrical projection}}, the vertical distances (along the same longitude line) are correct throughout. {space 8}{hline 27} {it:example do-file content} {hline 27} {cmd}{...} {* example_start - picard}{...} use "geo2xy_us_coor.dta", clear drop if inlist(_ID, 14, 39, 42) // Alaska, Puerto Rico, Hawaii geo2xy _Y _X , gen(lat lon) projection(picard, 38.889689) // the projected coordinates are the same as the original coordinates assert _Y == lat assert _X == lon // show the projection details and compute the plot's height return list local yheight = 6 * `r(aspect)' line lat lon, lwidth(vthin) lcolor(gray) 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(picard, replace) /// xline(-77.035279, lstyle(grid)) yline(38.889689, lstyle(grid)) {* example_end}{...} {txt}{...} {space 8}{hline 80} {space 8}{it:({stata geo2xy_run picard using geo2xy_picard.hlp, requires("geo2xy_us_coor.dta") preserve:click to run})} {* example_start - picard_t}{...} {* use "geo2xy_us_coor.dta", clear}{...} {* drop if inlist(_ID, 14, 39, 42) // Alaska, Puerto Rico, Hawaii}{...} {* }{...} {* geo2xy _Y _X , gen(lat lon) projection(picard, 38.889689) tissot}{...} {* }{...} {* // the projected coordinates are the same as the original coordinates}{...} {* assert _Y == lat}{...} {* assert _X == lon}{...} {* }{...} {* return list}{...} {* local yheight = 6 * `r(aspect)'}{...} {* }{...} {* line lat lon if !mi(_ID), lwidth(vthin) lcolor(gray) cmissing(n) ///}{...} {* || ///}{...} {* line lat lon 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) ///}{...} {* xline(-77.035279, lstyle(grid)) yline(38.889689, lstyle(grid))}{...} {* example_end}{...} {space 8}{it:({stata geo2xy_run picard_t using geo2xy_picard.hlp, requires("geo2xy_us_coor.dta") preserve:click to run with Tissot's indicatrices})} {title:Certification} {pstd} In the following example, the height of the map represents a 10 degree range of latitude (45-35). The width of the map represents 150 degrees of longitude (75 - -75) as measured at the mid-latitude point (40 degrees). {space 8}{hline 27} {it:example do-file content} {hline 27} {cmd}{...} {* example_start - certify}{...} clear input double(_Y _X) 35 -75 45 75 end // short-cut to convert decimal degrees to radians local d2r _pi / 180 local ar = (10 * `d2r') / (150 * `d2r' * cos(40 * `d2r')) geo2xy _Y _X, gen(y x) proj(picard) return list assert `ar' == `r(aspect)' assert _X == x assert _Y == y {* example_end}{...} {txt}{...} {space 8}{hline 80} {space 8}{it:({stata geo2xy_run certify using geo2xy_picard.hlp, preserve:click to run})} {title:References and further reading} {pstd} Snyder, John P., {it:Map projections: A working manual}. U.S. Geological Survey Professional Paper 1395, 1987. [{it:{browse "http://pubs.er.usgs.gov/publication/pp1395":download}}]. {pstd} Wikipedia, {it:{browse "http://en.wikipedia.org/wiki/Equirectangular_projection":Equirectangular projection}}. Web page, reviewed 1/31/2017. {title:Author} {pstd}Robert Picard{p_end} {pstd}picard@netbox.com{p_end}