{smcl} {* *! version 1.0.1 15aug2015}{...} {cmd:help geo2xy_proj} {hline} {title:Title} {phang} {bf:geo2xy} {hline 2} Supported map projections {pstd} The following projections are available: {synoptset 32 tabbed}{...} {synopthdr :Projection names} {synoptline} {synopt :{help geo2xy_proj##web_mercator:web_mercator}}Web Mercator projection - Spherical model{p_end} {synopt :{help geo2xy_proj##mercator_sphere:mercator_sphere}}Mercator projection - Spherical model{p_end} {synopt :{help geo2xy_proj##mercator:mercator}}Mercator projection - Ellipsoid model{p_end} {synopt :{help geo2xy_proj##equidistant_cylindrical:equidistant_cylindrical}}Equidistant Cylindrical projection - Spherical model{p_end} {synopt :{help geo2xy_proj##albers_sphere:albers_sphere}}Albers Equal-Area Conic Projection - Spherical model{p_end} {synopt :{help geo2xy_proj##albers:albers}}Albers Equal-Area Conic Projection - Ellipsoid model{p_end} {synopt :{help geo2xy_proj##picard:picard}}Picard's variant of the Equidistant Cylindrical projection - Spherical model{p_end} {synoptline} {p2colreset}{...} {marker web_mercator}{...} {title:Web Mercator projection - Spherical model} {pstd} {ul:{hi:Syntax}} {pstd} If no parameter is specified {p 8 16 2} {opt proj:ection}{hi:(web_mercator}) {pstd} If projection parameters are specified, all parameters must be specified and appear in the following order {p 8 16 2} {opt proj:ection}{hi:(web_mercator} {it:, zoom xtile ytile}) {synoptset 15 tabbed}{...} {synopthdr :Parameters} {synoptline} {synopt :{it:zoom}}zoom (default is 0){p_end} {synopt :{it:xtile}}name for new x tile id variable (default is no tile id generated){p_end} {synopt :{it:ytile}}name for new y tile id variable (default is no tile id generated){p_end} {synoptline} {p2colreset} {pstd} {ul:{hi:Description}} {pstd} The Web Mercator projection is used by most web mapping applications (Google Maps, Bing, etc.) and is {cmd:geo2xy}'s default if no projection is specified. This is a variant of a spherical Mercator projection. This implementation matches exactly what is done in Google Maps. A correctly scaled raster map of the projected coordinates can be overlaid on a Google base map and all points and features will align perfectly. {pstd} In a Mercator projection, longitudes and latitudes are projected as straight lines, perpendicular to each other. To compensate for the fact that longitudes are not parallel and converge at the pole, latitudes are stretched to match the local stretch in longitude. As longitudes converge towards the poles, the amount of stretch needed for latitudes to maintain proportions becomes infinite. {pstd} When Google implemented Google Maps, they decided that their maps would be rendered using tiles of 256x256 pixels each. A zoom level parameter decides the number of tiles needed to render the map of the whole world. At a zoom level of 0, the world fits in a single tile of 256x256 pixels. An increase of 1 in the zoom level doubles the number of tiles both horizontally and vertically. This imposes an aspect ratio of 1 on the map of the world. Since a Mercator projection stretches the y-coordinates towards infinity, Google's square map of the world truncates y-coordinates near the poles. On a unit sphere, the width of the projected equator for the whole world is 2 * _pi. That means that the projected latitudes need to be cut off at -_pi <= y <= _pi to stay within the bounds of a square tile. {pstd} Note that unless you are trying to overlay your maps with other maps, tile and pixel position are irrelevant and therefore the zoom level does not matter. In other words, if you plot the projected coordinates in Stata, the map produced will look exactly the same, no matter the zoom level. Furthermore, it will look exactly the same if you use a {help geo2xy_proj##mercator_sphere:spherical Mercator} projection instead. {pstd} With a Mercator projection, circles retain their shapes over small areas but their size increases towards the pole. {pstd} {ul:{hi:Spheroid and (x,y) coordinates units}} {pstd} This projection assumes that the geographic latitude and longitude describe locations on a sphere. The computations are performed on a unit sphere (radius of 1) and projected coordinates are returned in pixels from the map's origin (top-left) at the selected zoom level. {pstd} {ul:{hi:References and further reading}} {pstd} An overview of the Web Mercator projection can be found on the {browse "http://en.wikipedia.org/wiki/Web_Mercator":Web Mercator} entry of Wikipedia. {pstd} The following {browse "http://developers.google.com/maps/documentation/ios/tiles":Google Developers page} describes the basic structure of coordinates and tiles in Google Maps. {pstd} {ul:{hi:Example}} {pstd} This example is projected at a zoom level of 10. At that level, Google's map of the world is made up of 1024 * 1024 tiles (2^10) or 262,144 by 262,144 pixels. The top left of the map is (0,0) so {it:ylat} is negative. Of course {cmd:geo2xy} does not create a raster map, it just calculates the projected (x,y) coordinates in terms of pixel position and identifies the tile where the point is located. {cmd}. use "geo2xy_us_coor.dta", clear . keep if _ID == 44 . geo2xy _Y _X, gen(ylat xlon) project(web_mercator,10 xtile ytile) . return list . list in 1/10 . line ylat xlon, lwidth(vthin) lcolor(gray) cmissing(n) /// ylabel(minmax, nogrid) yscale(off) xlabel(minmax, nogrid) xscale(off) /// aspectratio(`r(aspect)') legend(off) . summarize{txt} {it:({stata geo2xy_examples ex_web_mercator:click to run})} {marker mercator_sphere}{...} {title:Mercator projection - Spherical model} {pstd} {ul:{hi:Syntax}} {p 8 16 2} {opt proj:ection}{hi:(mercator_sphere} [{it:, lon0}]) {synoptset 15 tabbed}{...} {synopthdr :Parameters} {synoptline} {synopt :{it:lon0}}projection's origin - longitude (default is mid-longitude){p_end} {synoptline} {p2colreset} {pstd} {ul:{hi:Description}} {pstd} This is exactly the same projection as the default {help geo2xy_proj##web_mercator:Web Mercator} projection but with no truncation of y-coordinates near the poles, no Google tiles, a different projection origin, and a different scaling of the (x,y) coordinates. {pstd} With a Mercator projection, circles retain their shapes over small areas but their size increases towards the pole. {pstd} {ul:{hi:Spheroid and (x,y) coordinates units}} {pstd} This projection assumes that the geographic latitude and longitude describe locations on a sphere. The computations are performed on a unit sphere (radius of 1) and projected coordinates are returned without further scaling. If you want the units to reflect distances on Earth, you need to multiply the coordinates by an appropriate radius (e.g. 6371 km, 3958.76 miles, 20,902,231 feet, etc.). {pstd} {ul:{hi:References and further reading}} {pstd} The equations for this projection come from page 41 of: {p 8 16 2} Snyder, John Parr. {it:Map projections--A working manual. No. 1395}. USGPO, 1987. {pstd} The numerical example at p. 266 was used to certify {bf:geo2xy}'s implementation. {pstd} An overview of the spherical Mercator projection can be found on the {browse "http://en.wikipedia.org/wiki/Mercator_projection":Mercator projection} entry of Wikipedia. {pstd} {ul:{hi:Example}} {pstd} This is a repeat of the {help geo2xy_proj##web_mercator:web_mercator} example above. The only difference is in the scaling of the (x,y) coordinates. Note the same aspect ratio. {cmd}. use "geo2xy_us_coor.dta", clear . keep if _ID == 44 . geo2xy _Y _X, gen(ylat xlon) project(mercator_sphere) tissot . return list . list in 1/10 . line ylat xlon if !mi(_ID), lwidth(vthin) lcolor(gray) cmissing(n) || /// line ylat xlon if mi(_ID), lwidth(vthin) lcolor(blue) cmissing(n) /// ylabel(minmax, nogrid) yscale(off) xlabel(minmax, nogrid) xscale(off) /// aspectratio(`r(aspect)') legend(off) . summarize{txt} {it:({stata geo2xy_examples ex_mercator_sphere:click to run})} {marker mercator}{...} {title:Mercator projection - Ellipsoid model} {pstd} {ul:{hi:Syntax}} {pstd} If no parameter is specified {p 8 16 2} {opt proj:ection}{hi:(mercator}) {pstd} If projection parameters are specified, all parameters must be specified and appear in the following order {p 8 16 2} {opt proj:ection}{hi:(mercator} {it:, a f lon0}) {synoptset 15 tabbed}{...} {synopthdr :Parameters} {synoptline} {synopt :{it:a}}semi-major axis of reference ellipsoid (default is 6378137){p_end} {synopt :{it:f}}inverse flattening of reference ellipsoid (default is 298.257223563){p_end} {synopt :{it:lon0}}projection's origin - longitude (default is mid-longitude){p_end} {synoptline} {p2colreset} {pstd} {ul:{hi:Description}} {pstd} With a Mercator projection, circles retain their shapes over small areas but their size increases towards the pole. {pstd} {ul:{hi:Spheroid and (x,y) coordinates units}} {pstd} This projection assumes that the geographic latitude and longitude describe locations on an ellipsoid. If parameters are not specified, {bf:geo2xy} uses the reference ellipsoid of the WGS 84 datum, the same that is used in GPS devices. The projected coordinates reflect the distance in meters from the projection's origin. {pstd} {ul:{hi:References and further reading}} {pstd} The equations for this projection come from page 44 of: {p 8 16 2} Snyder, John Parr. {it:Map projections--A working manual. No. 1395}. USGPO, 1987. {pstd} The numerical example at p. 267 was used to certify {bf:geo2xy}'s implementation. {pstd} {ul:{hi:Example}} {pstd} A finer point that is easy to overlook when working with geographic latitude and longitude is that the coordinates describe the location of a point on a reference spheroid. If you use a GPS to locate your position, the coordinates are on the reference ellipsoid of the WGS 84 datum. So if Ann Arbor's City Hall is located at 42.281716, -83.745176, then the distance between it and the equator along the same meridian is (requires {cmd:geodist} from SSC, {it:{stata ssc install geodist:click to install}}): . {stata geodist 42.281716 -83.745176 0 -83.745176} . {stata geodist 42.281716 -83.745176 0 -83.745176, sphere} {pstd} The difference in distance is amplified with a Mercator projection as the y-coordinates are increasingly stretched towards the poles. These differences matter to cartographers and so it's important to understand that a spherical Mercator map cannot be matched to one based on an ellipsoid. The spherical Mercator uses a unit sphere so the projected coordinates are converted to meters (the radius of WGS84 is 6378137 meters) to make them comparable to the ellipsoid version. The {help summarize} statements at the end of this example show that the difference in y-coordinates increases at higher latitudes: {cmd}. use "geo2xy_us_coor.dta", clear . keep if _ID == 44 . geo2xy _Y _X, gen(ys xs) project(mercator_sphere) . gen double xlons = xs * 6378137 . gen double ylats = ys * 6378137 . geo2xy _Y _X, gen(ylat xlon) project(mercator) . return list . list in 1/10 . line ylat xlon, lwidth(vthin) lcolor(gray) cmissing(n) || /// line ylats xlons, lwidth(vthin) lcolor(eltblue) cmissing(n) /// ylabel(minmax, nogrid) yscale(off) xlabel(minmax, nogrid) xscale(off) /// aspectratio(`r(aspect)') legend(off) . gen ydiff = ylats - ylat . summarize . summarize ydiff if _Y < 41.7 . summarize ydiff if _Y > 48.19{txt} {it:({stata geo2xy_examples ex_mercator:click to run})} {marker equidistant_cylindrical}{...} {title:Equidistant Cylindrical projection - Spherical model} {pstd} {ul:{hi:Syntax}} {pstd} If no parameter is specified {p 8 16 2} {opt proj:ection}{hi:(equidistant_cylindrical}) {pstd} If projection parameters are specified, all parameters must be specified and appear in the following order {p 8 16 2} {opt proj:ection}{hi:(equidistant_cylindrical} {it:, lat1 lon0}) {synoptset 15 tabbed}{...} {synopthdr :Parameters} {synoptline} {synopt :{it:lat1}}standard parallel (default is mid-latitude){p_end} {synopt :{it:lon0}}projection's origin (default is mid-longitude){p_end} {synoptline} {p2colreset} {pstd} {ul:{hi:Description}} {pstd} This is a very simple projection where latitudes are converted to radians (all vertical distances on the map are equal to the distance on a unit sphere) and longitudes are scaled such that the horizontal distance is accurate at the standard parallel ({it:lat1}). By default, this is mid-latitude in terms of all the points on the map. {pstd} At the top of the map, horizontal distances are exagerated while horizontal distances are compressed at the bottom of the map. {pstd} {ul:{hi:Spheroid and (x,y) coordinates units}} {pstd} This projection assumes that the geographic latitude and longitude describe locations on a sphere. The computations are performed on a unit sphere (radius of 1) and projected coordinates are returned without further scaling. If you want the units to reflect distances on Earth, you need to multiply the coordinates by an appropriate radius (e.g. 6371 km, 3958.76 miles, 20,902,231 feet, etc.). {pstd} {ul:{hi:References and further reading}} {pstd} The equations for this projection come from pages 90-91 of: {p 8 16 2} Snyder, John Parr. {it:Map projections--A working manual. No. 1395}. USGPO, 1987. {pstd} No certification on this one as Snyder puts it: "Numerical examples are omitted in the appendix, due to simplicity" (p. 91) {pstd} An overview of the Equidistant Cylindrical projection can be found on the {browse "http://en.wikipedia.org/wiki/Equirectangular_projection":Equirectangular projection} entry of Wikipedia. {pstd} {ul:{hi:Examples}} {pstd} This projection works well if the map area is not too large or near the poles: {cmd}. use "geo2xy_us_coor.dta", clear . keep if _ID == 44 . geo2xy _Y _X, gen(ylat xlon) project(equidistant_cylindrical) tissot . return list . line ylat xlon if !mi(_ID), lwidth(vthin) lcolor(gray) cmissing(n) || /// line ylat xlon if mi(_ID), lwidth(vthin) lcolor(blue) cmissing(n) /// ylabel(minmax, nogrid) yscale(off) xlabel(minmax, nogrid) xscale(off) /// aspectratio(`r(aspect)') legend(off){txt} {it:({stata geo2xy_examples ex_equi_mi:click to run})} {pstd} This is not the best choice for representing the state of Alaska (note that the longitudes must be shifted because the state spans the international date line): {cmd}. use "geo2xy_us_coor.dta", clear . keep if _ID == 37 . gen double _XX = cond(_X > 0, _X - 180, _X + 180) . geo2xy _Y _XX , gen(ylat xlon) proj(equidistant_cylindrical) tissot . line ylat xlon if !mi(_ID), lwidth(vthin) lcolor(gray) cmissing(n) || /// line ylat xlon if mi(_ID), lwidth(vthin) lcolor(eltblue) cmissing(n) /// ylabel(minmax, nogrid) yscale(off) xlabel(minmax, nogrid) xscale(off) /// aspectratio(`r(aspect)') legend(off){txt} {it:({stata geo2xy_examples ex_equi_ak:click to run})} {marker albers_sphere}{...} {title:Albers Equal-Area Conic Projection - Spherical model} {pstd} {ul:{hi:Syntax}} {pstd} If no parameter is specified {p 8 16 2} {opt proj:ection}{hi:(albers_sphere}) {pstd} If projection parameters are specified, all parameters must be specified and appear in the following order {p 8 16 2} {opt proj:ection}{hi:(albers_sphere} {it:, lat1 lat2 lat0 lon0}) {synoptset 15 tabbed}{...} {synopthdr :Parameters} {synoptline} {synopt :{it:lat1}}1st standard parallel (default is minlat + (maxlat - minlat) / 6){p_end} {synopt :{it:lat2}}2nd standard parallel (default is maxlat - (maxlat - minlat) / 6){p_end} {synopt :{it:lat0}}projection's origin (default is 0){p_end} {synopt :{it:lon0}}central meridian (default is mid-longitude){p_end} {synoptline} {p2colreset} {pstd} {ul:{hi:Description}} {pstd} This projection transfers locations on the spheroid to a cone that intersects the spheroid at the two standard parallels. The parallels are circles whose center is the intersection of the meridians. The algorithm to automatically select the first and second standard parallels follows Deetz and Adams's (1934, p. 91) suggestion (1/6 of meridional distance). {pstd} This is an equal area projection, which makes it particularly appropriate for choropleth maps. It is often used to represent the conterminous United States. {pstd} {ul:{hi:Spheroid and (x,y) coordinates units}} {pstd} This projection assumes that the geographic latitude and longitude describe locations on a sphere. The computations are performed on a unit sphere (radius of 1) and projected coordinates are returned without further scaling. If you want the units to reflect distances on Earth, you need to multiply the coordinates by an appropriate radius (e.g. 6371 km, 3958.76 miles, 20,902,231 feet, etc.). {pstd} {ul:{hi:References and further reading}} {pstd} Deetz, Charles Henry, and Oscar Sherman Adams. {it:Elements of map projection with applications to map and chart construction. No. 68}. US Government Printing Office, 1921. {pstd} The equations for this projection come from pages 100-101 of: {p 8 16 2} Snyder, John Parr. {it:Map projections--A working manual. No. 1395}. USGPO, 1987. {pstd} The numerical example at p. 291 was used to certify {bf:geo2xy}'s implementation. {pstd} An overview of the Albers Equal-Area Conic Projection can be found on the {browse "http://en.wikipedia.org/wiki/Albers_projection":Albers projection} entry of Wikipedia. {pstd} {ul:{hi:Example}} {pstd} {cmd: geo2xy} will automatically select appropriate values for the first and second parallel. For the conterminous United States, the standard parallels are usually set to 29.5 and 45.5. {cmd}. use "geo2xy_us_coor.dta", clear . drop if inlist(_ID, 55,54,32,28,29,37,40) // Alaska, Hawaii, territories . geo2xy _Y _X , gen(ylat xlon) proj(albers_sphere) . return list . sum _X . dis (r(min) + r(max)) / 2 . geo2xy _Y _X , gen(ylat2 xlon2) proj(albers_sphere, 29.5 45.5 0 -95.856482) . return list . line ylat2 xlon2, lwidth(vthin) lcolor(gray) cmissing(n) /// ylabel(minmax, nogrid) yscale(off) xlabel(minmax, nogrid) xscale(off) /// aspectratio(`r(aspect)') legend(off){txt} {it:({stata geo2xy_examples ex_albers_sphere:click to run})} {it:({stata geo2xy_examples ex_albers_sphere_t:click to run with Tissot's indicatrix})} {marker albers}{...} {title:Albers Equal-Area Conic Projection - Ellipsoid model} {pstd} {ul:{hi:Syntax}} {pstd} If no parameter is specified {p 8 16 2} {opt proj:ection}{hi:(albers}) {pstd} If projection parameters are specified, all parameters must be specified and appear in the following order {p 8 16 2} {opt proj:ection}{hi:(albers} {it:, a f lat1 lat2 lat0 lon0}) {synoptset 15 tabbed}{...} {synopthdr :Parameters} {synoptline} {synopt :{it:a}}semi-major axis of reference ellipsoid (default is 6378137){p_end} {synopt :{it:f}}inverse flattening of reference ellipsoid (default is 298.257223563){p_end} {synopt :{it:lat1}}1st standard parallel (default is minlat + (maxlat - minlat) / 6){p_end} {synopt :{it:lat2}}2nd standard parallel (default is maxlat - (maxlat - minlat) / 6){p_end} {synopt :{it:lat0}}projection's origin (default is 0){p_end} {synopt :{it:lon0}}central meridian (default is mid-longitude){p_end} {synoptline} {p2colreset} {pstd} {ul:{hi:Description}} {pstd} This projection transfers locations on the spheroid to a cone that intersects the spheroid at the two standard parallels. The parallels are circles whose center is the intersection of the meridians. The algorithm to automatically select the first and second standard parallels follows Deetz and Adams's (1934, p. 91) suggestion (1/6 of meridional distance). {pstd} This is an equal area projection, which makes it particularly appropriate for choropleth maps. It is often used to represent the conterminous United States. {pstd} {ul:{hi:Spheroid and (x,y) coordinates units}} {pstd} This projection assumes that the geographic latitude and longitude describe locations on an ellipsoid. If parameters are not specified, {bf:geo2xy} uses the reference ellipsoid of the WGS 84 datum, the same that is used in GPS devices. The projected coordinates reflect distances in meters. {pstd} {ul:{hi:References and further reading}} {pstd} Deetz, Charles Henry, and Oscar Sherman Adams. {it:Elements of map projection with applications to map and chart construction. No. 68}. US Government Printing Office, 1921. {pstd} The equations for this projection come from pages 101-102 of: {p 8 16 2} Snyder, John Parr. {it:Map projections--A working manual. No. 1395}. USGPO, 1987. {pstd} The numerical example at p. 292-293 was used to certify {bf:geo2xy}'s implementation. {pstd} An overview of the Albers Equal-Area Conic Projection can be found on the {browse "http://en.wikipedia.org/wiki/Albers_projection":Albers projection} entry of Wikipedia. {pstd} {ul:{hi:Example}} {pstd} {cmd: geo2xy} will automatically select appropriate values for the first and second parallel. For the conterminous United States, the standard parallels are usually set to 29.5 and 45.5. {cmd}. use "geo2xy_us_coor.dta", clear . drop if inlist(_ID, 55,54,32,28,29,37,40) // Alaska, Hawaii, territories . geo2xy _Y _X , gen(ylat xlon) proj(albers) . return list . sum _X . dis (r(min) + r(max)) / 2 . geo2xy _Y _X , gen(ylat2 xlon2) proj(albers, 6378137 298.257223563 29.5 45.5 0 -95.856482) . return list . line ylat2 xlon2, lwidth(vthin) lcolor(gray) cmissing(n) /// ylabel(minmax, nogrid) yscale(off) xlabel(minmax, nogrid) xscale(off) /// aspectratio(`r(aspect)') legend(off){txt} {it:({stata geo2xy_examples ex_albers:click to run})} {it:({stata geo2xy_examples ex_albers_t:click to run with Tissot's indicatrix})} {marker picard}{...} {title:Picard's Equidistant Cylindrical projection} {pstd} {ul:{hi:Syntax}} {p 8 16 2} {opt proj:ection}{hi:(picard} [{it:, lat1}]) {synoptset 15 tabbed}{...} {synopthdr :Parameters} {synoptline} {synopt :{it:lat1}}standard parallel (default is mid-latitude){p_end} {synoptline} {p2colreset} {pstd} {ul:{hi:Description}} {pstd} On a unit sphere, the distance as measured on 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 {help geo2xy_proj##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. {pstd} {ul:{hi: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. {pstd} {ul:{hi:Example}} {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 {help geo2xy_proj##equidistant_cylindrical:Equidistant Cylindrical projection}, the vertical distances (along the same longitude line) are correct throughout. {cmd}. use "geo2xy_us_coor.dta", clear . drop if inlist(_ID, 55,54,32,28,29,37,40) // Alaska, Hawaii, territories . geo2xy _Y _X , gen(lat lon) proj(picard, 38.889689) . return list . line lat lon, lwidth(vthin) lcolor(gray) cmissing(n) /// ylabel(minmax, nogrid) yscale(off) xlabel(minmax, nogrid) xscale(off) /// aspectratio(`r(aspect)') legend(off) /// xline(-77.035279, lstyle(grid)) yline(38.889689, lstyle(grid)){txt} {it:({stata geo2xy_examples ex_picard:click to run})} {it:({stata geo2xy_examples ex_picard_t:click to run with Tissot's indicatrix})}