{smcl}
{* *! version 1.0.0  31jan2017}{...}
{cmd:help geo2xy_web_mercator}
{hline}

{title:Title}

{phang}
{cmd:geo2xy} {hline 2} Convert latitude and longitude to cartesian (x,y) coordinates


{title:Map projection}

{phang}
Web mercator projection - spherical model (Google Maps, Bing, etc.)


{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 web_mercator}
[,{help geo2xy##proj_name:proj_opts}]{cmd:)}
{opt ti:ssot} ]


{synoptset 40 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(web_mercator [,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: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}
Projection parameters are optional. To specify some, 
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})


{title:Description}

{pstd}
The Web Mercator projection is used by most
web mapping applications (Google Maps, Bing, etc.).
It is the same as a standard Mercator projection
except that it uses a sphere as a model for the earth
instead of an ellipsoid. 

{pstd}
This projection, developed by Google for their Google Maps service,
returns projected
coordinates in map pixel coordinates and optionally with map tile coordinates.

{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 poles, 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 truncates
the map of the world near the poles. Points at latitudes beyond 85.051129 degree 
north or south have missing y-coordinates.

{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.


{title: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.


{title:Examples}

{pstd}
These examples require {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}
Since this is {cmd:geo2xy}'s default projection, there's no need
to specify a projection. 

{space 8}{hline 27} {it:example do-file content} {hline 27}
{cmd}{...}
{* example_start - web_mercator_us}{...}
	use "geo2xy_us_coor.dta", clear
	drop if inlist(_ID, 14, 39, 42) // Alaska, Puerto Rico, Hawaii

	geo2xy _Y _X, gen(ylat xlon)
	
	// show the projection details and compute the plot's height
	return list
	local yheight = 6 * `r(aspect)'
	
	line ylat xlon, 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(web_mercator_us, replace)
			
	summarize
{* example_end}{...}
{txt}{...}
{space 8}{hline 80}
{space 8}{it:({stata geo2xy_run web_mercator_us using geo2xy_web_mercator.hlp, requires("geo2xy_us_coor.dta") preserve:click to run})}

{pstd}
Redo the same map, this time 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.

{space 8}{hline 27} {it:example do-file content} {hline 27}
{cmd}{...}
{* example_start - web_mercator_us_z10}{...}
	use "geo2xy_us_coor.dta", clear
	drop if inlist(_ID, 14, 39, 42) // Alaska, Puerto Rico, Hawaii
	
	geo2xy _Y _X, gen(ylat xlon) projection(web_mercator, 10 xtile ytile)
	
	// show the projection details and compute the plot's height
	return list
	local yheight = 6 * `r(aspect)'
	
	line ylat xlon, 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(web_mercator_us_z10, replace)
			
	summarize
{* example_end}{...}
{txt}{...}
{space 8}{hline 80}
{space 8}{it:({stata geo2xy_run web_mercator_us_z10 using geo2xy_web_mercator.hlp, requires("geo2xy_us_coor.dta") preserve:click to run})}

{pstd}
As you can see, the projected coordinates are different but the 
maps produced by the two examples above are identical.

{pstd}
Use the {opt tissot} option to illustrate the distortions created by
this projection:

{space 8}{hline 27} {it:example do-file content} {hline 27}
{cmd}{...}
{* example_start - web_mercator_tissot}{...}
	use "geo2xy_us_coor.dta", clear
	drop if inlist(_ID, 14, 39, 42) // Alaska, Puerto Rico, Hawaii

	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(web_mercator_tissot, replace)
{* example_end}{...}
{txt}{...}
{space 8}{hline 80}
{space 8}{it:({stata geo2xy_run web_mercator_tissot using geo2xy_web_mercator.hlp, requires("geo2xy_us_coor.dta") preserve:click to run})}


{title:Certification}

{pstd}
Google Maps offers a 
{it:{browse "https://developers.google.com/maps/documentation/javascript/examples/map-coordinates":JavaScript utility}}
that shows both the pixel
and tile coordinates of a location in Google Maps.
Their current example uses a point in Chicago, IL at latitude 41.85
and longitude -87.650. If you zoom in or out, the caption
shows the pixel and tile coordinates. For example:

{pmore}
Zoom level: 14  Pixel Coordinate: (1075955, 1559345)  Tile Coordinate: (4202, 6091)

{pstd}
The following code confirms that {cmd:geo2xy} generates the same projected
pixel and tile coordinates.
Note that since
the map origin (0,0) is the top-left corner, y-coordinates are negative.

{space 8}{hline 27} {it:example do-file content} {hline 27}
{cmd}{...}
{* example_start - certify}{...}
	clear
	input double(_Y _X)
	41.85 -87.650
	end

	geo2xy _Y _X, gen(y x) projection(web_mercator, 14 xtile ytile)
	return list
	list
	
	assert int(x) == 1075955
	assert int(y) == -1559345
	assert xtile == 4202
	assert ytile == 6091
{* example_end}{...}
{txt}{...}
{space 8}{hline 80}
{space 8}{it:({stata geo2xy_run certify using geo2xy_web_mercator.hlp, preserve:click to run})}


{title:References and further reading}

{pstd}
Google, {it:{browse "http://developers.google.com/maps/documentation/ios/tiles":Tile Layers}}.
Web page, reviewed 1/31/2017.

{pstd}
Microsoft, {it:{browse "https://msdn.microsoft.com/en-us/library/bb259689.aspx":Bing Maps Tile System}}.
Web page, reviewed 1/31/2017.

{pstd}
Wikipedia, {it:{browse "http://en.wikipedia.org/wiki/Web_Mercator":Web Mercator}}.
Web page, reviewed 1/31/2017.

{pstd}
Since this is a simple variant of a Mercator with a spherical model, see {cmd:geo2xy}'s
{help geo2xy_mercator_sphere:Mercator projection - spherical model} 
for more details about the underlying projection.


{title:Author}

{pstd}Robert Picard{p_end}
{pstd}picard@netbox.com{p_end}