*! version 1.0.0 24dec2016 Michael Stepner program define geo2xy_lambert_sphere, rclass /* ------------------------------------------------------------------------------- Lambert Conformal 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 106-107. Certified using numerical example at pp. 295 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. 107, equation 15-3 tempname n scalar `n' = ln(cos(`phi1')/cos(`phi2')) / ln( tan(_pi/4 + `phi2'/2) / tan(_pi/4 + `phi1'/2) ) // Snyder, p. 107, equation 15-2 tempname F scalar `F' = cos(`phi1') * tan(_pi/4 + `phi1'/2)^`n' / `n' // Snyder, p. 106, equation 15-1a tempname rho0 scalar `rho0' = `F' / (tan(_pi/4 + `phi0'/2)^`n') // Snyder, p. 106, equation 14-4 tempvar theta gen double `theta' = `n' * (`lon' * `d2r' - `lambda0') // Snyder, p. 106, equation 15-1 tempvar rho gen double `rho' = `F' / (tan(_pi/4 + (`lat' * `d2r')/2)^`n') local genrep = cond("`lon'" == "`x'","replace","gen double") nobreak { // Snyder, p. 106, equation 14-1 `genrep' `x' = `rho' * sin(`theta') if `touse' // Snyder, p. 106, 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 = "Lambert Conformal 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