*! version 1.0.1  20jan2017 Robert Picard, picard@netbox.com
program define geo2xy, rclass

	version 9.2

	syntax varlist(min=2 max=2 numeric) [if] [in], 	///
		[ ///
		GENerate(namelist min=2 max=2) ///
		PROJection(string) ///
		TIssot ///
		replace ///
		]

	marksample touse
	qui count if `touse'
	if r(N) == 0 error 2000

	tokenize `varlist'
	local lat `1'
	local lon `2'

	sum `lat' if `touse', meanonly
	if r(max) > 90 | r(min) < -90 {
		dis as err "latitude `lat' must be between -90 and 90"
		exit 198
	}

	sum `lon' if `touse', meanonly
	if r(max) > 180 | r(min) < -180 {
		dis as err "longitude `lon' must be between -180 and 180"
		exit 198
	}

	if "`generate'" == "" & "`replace'" == "" {
		dis as err "you must either specify the -replace- option or " ///
			"use the -generate()- option for the x,y coordinates variables"
		exit 198
	}

	if "`generate'" == "" {
		local yvar `lat'
		local xvar `lon'
	}
	else {
		tokenize `generate'
		local yvar `1'
		local xvar `2'
	}

	if "`tissot'" != "" geo2xy_tissot `touse' `lat' `lon'

	if "`projection'" == "" local projection "web_mercator"

	gettoken which_proj proj_options : projection, parse(",")
	gettoken comma proj_options : proj_options, parse(",")

	geo2xy_`which_proj' `touse' `lat' `lon' `yvar' `xvar' `proj_options'

	return add

end


program geo2xy_tissot

	args touse lat lon

	tempname d2r r2d
	scalar `d2r' = _pi / 180
	scalar `r2d' = 180 / _pi

	// get coordinates of map bounds
	sum `lat' if `touse', meanonly
	tempname latmin latmax
	scalar `latmin' = r(min) * `d2r'
	scalar `latmax' = r(max) * `d2r'

	sum `lon' if `touse', meanonly
	tempname lonmin lonmax ratio
	scalar `lonmin' = r(min) * `d2r'
	scalar `lonmax' = r(max) * `d2r'


	// get the map ratio at mid-latitude
	scalar `ratio' = (`latmax' - `latmin') / ///
		(cos((`latmax' + `latmin') / 2) * (`lonmax' - `lonmin'))


	// decide on a number of indicatrix, use odd numbers to align
	// with mid-latitude and mid-longitude
	local nlat = 5
	local nlon = 2 * floor(`nlat' / `ratio' / 2) + 1

	preserve

	qui {

	drop _all

		// coordinates of indicatrix
		tempvar ylat xlon
		set obs `nlat'
		gen double `ylat' = `latmin' + _n * (`latmax' - `latmin') / (`nlat' + 1)
		expand `nlon'
		bys `ylat': gen double `xlon' = `lonmin' + _n * (`lonmax' - `lonmin') / (`nlon' + 1)

		// the radius to use is in proportion to the number of points
		local r = cos((`latmax' + `latmin') / 2) * (`lonmax' - `lonmin') / (`nlon' + 1) / 3

		// create a polygon of 200 points; first is missing, last loops back to first
		expand 202

		tempvar bearing
		bysort `ylat' `xlon': gen double `bearing' = 2 * _pi * (_n - 2) / (_N-2) if _n > 1

		/*
		The equation for calculating a destination point using an initial bearing and
		distance comes from www.movable-type.co.uk/scripts/latlong.html by Chris Veness
		*/
		gen double `lat' = asin(sin(`ylat') * cos(`r') + ///
							cos(`ylat') * sin(`r') * cos(`bearing'))

		gen double `lon' = `xlon' + atan2(sin(`bearing') * sin(`r') * cos(`ylat'), ///
											cos(`r') - sin(`ylat') * sin(`lat'))

		replace `lat' = `lat' * `r2d'
		replace `lon' = `lon' * `r2d'

		replace `lon' = `lon' - 360 if `lon' > 180
		replace `lon' = `lon' + 360 if `lon' < -180

		keep `lat' `lon'
		tempfile hold
		qui save "`hold'"

		restore

		append using "`hold'"

	}

end


program define geo2xy_web_mercator, rclass

/*
-------------------------------------------------------------------------------

Web Mercator
============

Google Maps and other web mapping application use a spherical mercator
projection (http://en.wikipedia.org/wiki/Web_Mercator).

This implementation follows Google's implementation (haven't checked what
the others do), where the projected map's origin (0,0) is the top-left
corner and the projected latitudes are truncated so that whole world fits,
at a zoom level of 0, in a 256x256 tile, i.e. an aspect ratio of 1.
(https://developers.google.com/maps/documentation/ios/tiles)

By default, the map coordinates are returned at zoom level 0 (zoom level
does not matter unless you are trying to overlay tiles with Google's tiles).

If a zoom level is specified, the xtile and ytile identifier are also
returned.


Projection arguments:

  zoom  = zoom level
  xtile = x tile id
  ytile = y tile id


Defaults if no projection arguments supplied

  zoom  = 0
  xtile =
  ytile =

-------------------------------------------------------------------------------
*/

	args touse lat lon y x zoom xtile ytile whatsthis

	if !mi("`whatsthis'") {
		dis as err "expected arguments are: lon0"
		exit 198
	}

	// default zoom level
	if "`zoom'" == ""  local zoom 0

	// at zoom level 0, the world fits in a single 256x256 pixel tile
	// an increase of 1 for zoom level doubles x and doubles y
	tempname zoomfactor
	scalar `zoomfactor' = 256 / (2 * _pi) * 2^`zoom'

	// offset that will shift geodetic longitudes so that -180 == 0
	local lon0 -180

	// start with a regular spherical Mercator projection
	geo2xy_mercator_sphere `touse' `lat' `lon' `y' `x' `lon0'

	qui {

		// truncate to maintain Google's 1 to 1 aspect ratior at zoom level 0
		// where x == 2 * _pi so -_pi <= y <= _pi
		replace `y' = . if abs(`y') > _pi

		// shift the projected latitude so that top-left == 0,0 and scale
		replace `y' = `zoomfactor' * (`y' - _pi)  if `touse'

		// the longitudes have already been shifted
		replace `x' = `zoomfactor' * `x'  if `touse'

		// if arguments are specified, generate tile identifiers
		if "`xtile'" != "" {

			gen long `xtile' = int(`x'/256)
			gen long `ytile' = int(-`y'/256)

		}

	}

	return local zoom `zoom'
	return local pname = "Web Mercator"
	return local model = "Unit Sphere"

	sum `y' if `touse', meanonly
	tempname height
	scalar `height' = r(max) - r(min)
	sum `x' if `touse', meanonly
	return local aspect = `height' / (r(max) - r(min))

end


program define geo2xy_mercator_sphere, rclass

/*
-------------------------------------------------------------------------------

Mercator (sphere)
=================

Source:

  Snyder, John Parr. Map projections--A working manual. No. 1395. USGPO, 1987.

Available from http://pubs.usgs.gov/pp/1395/report.pdf.

Formulas on pages 41. Certified using numerical example at pp. 266

Projection arguments:

  lon0  = projection's origin

Defaults if no projection arguments supplied

  lon0  = relative; set to mid longitude range

-------------------------------------------------------------------------------
*/

	args touse lat lon y x lon0 whatsthis

	if !mi("`whatsthis'") {
		dis as err "expected arguments are: lon0"
		exit 198
	}

	if "`lon0'" == "" {

		sum `lon' if `touse', meanonly
		local lon0 =  (r(min) + r(max)) / 2

	}

	tempname d2r
	scalar `d2r' = _pi / 180

	tempname lambda0
	scalar `lambda0' = `lon0' * `d2r'

	qui nobreak {

		local genrep = cond("`lon'" == "`x'","replace","gen double")

		// Snyder, p. 41, equation 7-1; R = 1
		`genrep' `x' = (`lon' * `d2r' - `lambda0') if `touse'

		// Snyder, p. 41, equation 7-2; R = 1
		`genrep' `y' = log(tan(_pi / 4 + `lat' * `d2r' / 2)) if `touse'

	}

	return local lon0 `lon0'
	return local pname = "Mercator"
	return local model = "Unit Sphere"

	sum `y' if `touse', meanonly
	tempname height
	scalar `height' = r(max) - r(min)
	sum `x' if `touse', meanonly
	return local aspect = `height' / (r(max) - r(min))

end


program define geo2xy_mercator, rclass

/*
-------------------------------------------------------------------------------

Mercator (ellipsoid)
====================

Source:

  Snyder, John Parr. Map projections--A working manual. No. 1395. USGPO, 1987.

Available from http://pubs.usgs.gov/pp/1395/report.pdf.

Formulas on pages 44. Certified using numerical example at pp. 267

Projection arguments:

a     = semi-major axis of reference ellipsoid
f     = inverse flattening of reference ellipsoid
lon0  = projection's origin

Defaults if no projection arguments supplied

a     = 6378137 semi-major axis of WGS84
f     = 298.257223563 inverse flattening of WGS84
lon0  = relative; set to mid longitude range

-------------------------------------------------------------------------------
*/

	args touse lat lon y x a f lon0 whatsthis

	if !mi("`whatsthis'") | (!mi("`a'") & mi("`lon0'")) {
		dis as err "expected arguments are: a f lon0"
		exit 198
	}

	if "`a'" == "" {

		local a 6378137
		local f 298.257223563

		sum `lon' if `touse', meanonly
		local lon0 =  (r(min) + r(max)) / 2

	}

	// semi-minor axis
	tempname b
	scalar `b' = `a' - `a' / `f'

	// eccentricity, see Snyder, p. 13
	tempname e e2
	scalar `e2' = 2 * (1/`f') - (1/`f')^2
	scalar `e' = sqrt(`e2')

	tempname d2r
	scalar `d2r' = _pi / 180

	tempname lambda0
	scalar `lambda0' = `lon0' * `d2r'

	qui nobreak {

		local genrep = cond("`lon'" == "`x'","replace","gen double")

		// Snyder, p. 44, equation 7-6
		`genrep' `x' = `a' * (`lon' * `d2r' - `lambda0')  if `touse'

		// Snyder, p. 44, equation 7-7
		`genrep' `y' = `a' * log( ///
			tan(_pi / 4 + `lat' * `d2r' / 2) * ///
			((1 -`e' * sin(`lat' * `d2r')) / ///
			(1 + `e' * sin(`lat' * `d2r')))^(`e' / 2)) if `touse'

	}

	return local a `a'
	return local f `f'
	return local lon0 `lon0'
	return local pname = "Mercator"
	return local model = "Ellipsoid (`a',`f')"

	sum `y' if `touse', meanonly
	tempname height
	scalar `height' = r(max) - r(min)
	sum `x' if `touse', meanonly
	return local aspect = `height' / (r(max) - r(min))

end


program define geo2xy_equidistant_cylindrical, rclass

/*
-------------------------------------------------------------------------------

Equidistant Cylindrical
=======================

Source:

  Snyder, John Parr. Map projections--A working manual. No. 1395. USGPO, 1987.

Available from http://pubs.usgs.gov/pp/1395/report.pdf.

Formulas on pages 90-91.

No certification as "Numerical examples are omitted in the appendix, due to
simplicity" (p. 91)

Projection arguments:

lat1  = standard parallel
lon0  = projection's origin

Defaults if no projection arguments supplied

lat1  = relative; mid-latitude
lon0  = relative; set to mid-longitude range

-------------------------------------------------------------------------------
*/

	args touse lat lon y x lat1 lon0 whatsthis

	if !mi("`whatsthis'") | (!mi("`lat1'") & mi("`lon0'")) {
		dis as err "expected arguments are: lat1 lon0"
		exit 198
	}

	// if no standard parallel is specified, use mid-latitude
	if "`lat1'" == "" {

		sum `lat' if `touse', meanonly
		local lat1 = (r(min) + r(max)) / 2

		sum `lon' if `touse', meanonly
		local lon0 =  (r(min) + r(max)) / 2

	}

	tempname d2r
	scalar `d2r' = _pi / 180

	tempname phi1 lambda0
	scalar `lambda0' = `lon0' * `d2r'
	scalar `phi1' = `lat1' * `d2r'

	qui nobreak {

		local genrep = cond("`lon'" == "`x'","replace","gen double")

		// Snyder, p. 91, equation 12-1; R = 1
		`genrep' `x' = (`lon' * `d2r' - `lambda0') * cos(`phi1') if `touse'

		// Snyder, p. 91, equation 12-2; R = 1
		`genrep' `y' = `lat' * `d2r' if `touse'


	}

	return local lat1 `lat1'
	return local lon0 `lon0'
	return local pname = "Equidistant Cylindrical"
	return local model = "Unit Sphere"

	sum `y' if `touse', meanonly
	tempname height
	scalar `height' = r(max) - r(min)
	sum `x' if `touse', meanonly
	return local aspect = `height' / (r(max) - r(min))

end


program define geo2xy_albers_sphere, rclass

/*
-------------------------------------------------------------------------------

Albers Equal-Area Conic Projection (sphere)
===========================================

Source:

  Snyder, John Parr. Map projections--A working manual. No. 1395. USGPO, 1987.

Available from http://pubs.usgs.gov/pp/1395/report.pdf.

Formulas on pages 100-101.

Certified using numerical example at pp. 291

Projection arguments:

lat1  = 1st standard parallel
lat2  = 2nd standard parallel
lat0  = projection's origin
lon0  = central meridian

Defaults if no projection arguments supplied

lat1  = relative; minlat + (maxlat - minlat) / 6
lat2  = relative; maxlat - (maxlat - minlat) / 6
lat0  = 0
lon0  = relative; set to mid-longitude range

-------------------------------------------------------------------------------
*/

	args touse lat lon y x lat1 lat2 lat0 lon0 whatsthis

	if !mi("`whatsthis'") | (!mi("`lat1'") & mi("`lon0'")) {
		dis as err "expected arguments are: lat1 lat2 lat0 lon0"
		exit 198
	}

	if "`lat1'" == "" {

		// follow Deetz & Adams 1/6 suggestion, Snyder, p. 99
		sum `lat' if `touse', meanonly
		local lat1 = r(min) + (r(max) - r(min)) / 6
		local lat2 = r(max) - (r(max) - r(min)) / 6

		sum `lon' if `touse', meanonly
		local lon0 =  (r(min) + r(max)) / 2

		local lat0 0

	}


	tempname d2r
	scalar `d2r' = _pi / 180

	tempname phi0 phi1 phi2 lambda0
	scalar `phi0' = `lat0' * `d2r'
	scalar `lambda0' = `lon0' * `d2r'
	scalar `phi1' = `lat1' * `d2r'
	scalar `phi2' = `lat2' * `d2r'

	qui {

		// Snyder, p. 100, equation 14-6
		tempname n
		scalar `n' = (sin(`phi1') + sin(`phi2')) / 2

		// Snyder, p. 100, equation 14-5
		tempname C
		scalar `C' = cos(`phi1')^2 + 2 * `n' * sin(`phi1')

		// Snyder, p. 100, equation 14-3a
		tempname rho0
		scalar `rho0' = sqrt(`C' - 2 * `n' * sin(`phi0')) / `n'

		// Snyder, p. 100, equation 14-4
		tempvar theta
		gen double `theta' = `n' * (`lon' * `d2r' - `lambda0')

		// Snyder, p. 100, equation 14-3
		tempvar rho
		gen double `rho' = sqrt(`C' - 2 * `n' * sin(`lat' * `d2r')) / `n'

		local genrep = cond("`lon'" == "`x'","replace","gen double")

		nobreak {

			// Snyder, p. 100, equation 14-1
			`genrep' `x' = `rho' * sin(`theta') if `touse'

			// Snyder, p. 100, equation 14-2
			`genrep' `y' = `rho0' - `rho' * cos(`theta') if `touse'


		}
	}

	return local lat0 `lat0'
	return local lat1 `lat1'
	return local lat2 `lat2'
	return local lon0 `lon0'
	return local pname = "Albers Equal-Area Conic Projection"
	return local model = "Unit Sphere"

	sum `y' if `touse', meanonly
	tempname height
	scalar `height' = r(max) - r(min)
	sum `x' if `touse', meanonly
	return local aspect = `height' / (r(max) - r(min))

end


program define geo2xy_albers, rclass

/*
-------------------------------------------------------------------------------

Albers Equal-Area Conic Projection (ellipsoid)
==============================================

Source:

  Snyder, John Parr. Map projections--A working manual. No. 1395. USGPO, 1987.

Available from http://pubs.usgs.gov/pp/1395/report.pdf.

Formulas on pages 101-102.

Certified using numerical example at pp. 292-293

Projection arguments:

a     = semi-major axis of reference ellipsoid
f     = inverse flattening of reference ellipsoid
lat1  = 1st standard parallel
lat2  = 2nd standard parallel
lat0  = projection's origin
lon0  = central meridian

Defaults if no projection arguments supplied

a     = 6378137 semi-major axis of WGS84
f     = 298.257223563 inverse flattening of WGS84
lat1  = relative; minlat + (maxlat - minlat) / 6
lat2  = relative; maxlat - (maxlat - minlat) / 6
lat0  = 0
lon0  = relative; set to mid-longitude range

-------------------------------------------------------------------------------
*/

	args touse lat lon y x a f lat1 lat2 lat0 lon0 whatsthis

	if !mi("`whatsthis'") | (!mi("`a'") & mi("`lon0'")) {
		dis as err "expected arguments are: a f lat1 lat2 lat0 lon0"
		exit 198
	}

	if "`a'" == "" {
		local a 6378137
		local f 298.257223563

		// follow Deetz & Adams 1/6 suggestion, Snyder, p. 99
		sum `lat' if `touse', meanonly
		local lat1 = r(min) + (r(max) - r(min)) / 6
		local lat2 = r(max) - (r(max) - r(min)) / 6

		sum `lon' if `touse', meanonly
		local lon0 =  (r(min) + r(max)) / 2

		local lat0 0
	}

	// semi-minor axis
	tempname b
	scalar `b' = `a' - `a' / `f'

	// eccentricity, see Snyder, p. 13
	tempname e e2
	scalar `e2' = 2 * (1/`f') - (1/`f')^2
	scalar `e' = sqrt(`e2')

	tempname d2r
	scalar `d2r' = _pi / 180

	tempname phi0 phi1 phi2 lambda0
	scalar `phi0' = `lat0' * `d2r'
	scalar `phi1' = `lat1' * `d2r'
	scalar `phi2' = `lat2' * `d2r'
	scalar `lambda0' = `lon0' * `d2r'

	qui {
		tempvar phi
		gen double `phi' = `lat' * `d2r'

		// Snyder, p. 101, equation 3-12 with subscripts 0, 1, 2, and none
		foreach i in 0 1 2 {
			tempname q`i'
			scalar `q`i'' = (1 - `e2') * ///
				(sin(`phi`i'') / (1 - `e2' * sin(`phi`i'')^2) - ///
				log((1 - `e' * sin(`phi`i'')) / (1 + `e' * sin(`phi`i''))) / (2 * `e'))
		}
		tempvar q
		gen double `q' = (1 - `e2') * ///
			(sin(`phi') / (1 - `e2' * sin(`phi')^2) - ///
			log((1 - `e' * sin(`phi')) / (1 + `e' * sin(`phi'))) / (2 * `e'))

		// Snyder, p. 101, equation 14-15 with subscripts 1 and 2
		tempname m1 m2
		scalar `m1' = cos(`phi1') / sqrt(1 - `e2' * sin(`phi1')^2)
		scalar `m2' = cos(`phi2') / sqrt(1 - `e2' * sin(`phi2')^2)

		// Snyder, p. 101, equation 14-14
		tempname n
		scalar `n' = (`m1'^2 - `m2'^2) / (`q2' - `q1')

		// Snyder, p. 101, equation 14-13
		tempname C
		scalar `C' = `m1'^2 + `n' * `q1'

		// Snyder, p. 101, equation 14-12a
		tempname rho0
		scalar `rho0' = `a' * sqrt(`C' - `n' * `q0') / `n'

		// Snyder, p. 101, equation 14-4
		tempvar theta
		gen double `theta' = `n' * (`lon' * `d2r' - `lambda0')

		// Snyder, p. 101, equation 14-12
		tempvar rho
		gen double `rho' = `a' * sqrt(`C' - `n' * `q') / `n'

		local genrep = cond("`lon'" == "`x'","replace","gen double")

		nobreak {

			// Snyder, p. 101, equation 14-1
			`genrep' `x' =  `rho' * sin(`theta') if `touse'

			// Snyder, p. 101, equation 14-2
			`genrep' `y' = `rho0' - `rho' * cos(`theta') if `touse'


		}
	}

	return local a `a'
	return local f `f'
	return local lat0 `lat0'
	return local lat1 `lat1'
	return local lat2 `lat2'
	return local lon0 `lon0'
	return local pname = "Albers Equal-Area Conic Projection"
	return local model = "Ellipsoid (`a',`f')"

	sum `y' if `touse', meanonly
	tempname height
	scalar `height' = r(max) - r(min)
	sum `x' if `touse', meanonly
	return local aspect = `height' / (r(max) - r(min))

end


program define geo2xy_picard, rclass

/*
-------------------------------------------------------------------------------

picard
======

The geographic coordinates remain unchanged. The aspect ratio is computed
using the distance at the standard parallel between min and max longitude.
When plotting the geographic coordinates using this aspect ratio will
produce exactly the same map as would be produced if the Equidistant
Cylindrical was used. The advantage is that you can add plot features
based on coordinates (lines, coordinates, etc.).

Projection arguments:

lat1  = standard parallel

Defaults if no projection arguments supplied

lat1  = relative; mid-latitude

-------------------------------------------------------------------------------
*/

	args touse lat lon y x lat1 whatsthis

	if !mi("`whatsthis'") {
		dis as err "expected arguments are: lat1"
		exit 198
	}

	if "`lat'" != "`y'" {
		qui clonevar `y' = `lat'
		qui clonevar `x' = `lon'
	}
	else {
		dis as txt "note: this projection does not change the coordinates"
	}

	// if no standard parallel is specified, use mid-latitude
	if "`lat1'" == "" {

		sum `lat' if `touse', meanonly
		local lat1 = (r(min) + r(max)) / 2

	}

	tempname d2r
	scalar `d2r' = _pi / 180

	return local lat1 `lat1'
	return local pname = "Picard"
	return local model = "Geographic Latitude & Longitude"

	sum `lat' if `touse', meanonly
	tempname height
	scalar `height' = (r(max) - r(min)) * `d2r'

	sum `lon' if `touse', meanonly
	return local aspect = `height' / (cos(`lat1' * `d2r') * (r(max) * `d2r' - r(min) * `d2r'))

end