Projections with geo2xy

[John MacDonald]

I'm having a strange problem with geo2xy that I can't figure out. I have a set of coordinates and a shape file that are both in WGS1984. When I used SPMAP they look fine. I can see the x y locations on top of the polygon of Philadelphia block groups.

here is the basic code example:
use "C:~phillyblockgsdb.dta", clear

spmap using "C:~phillyblockgscoord.dta", id(id) ///
fcolor(eggshell) ocolor(black) osize(thin) ///
point(data("C:~cluster_locations.dta") x(longitude_x) y (latitude_y) ///
size(*0.5) fcolor(green blue black) ocolor(white) ///
by (arm) ///
osize(vvthin)) ///
title ("Vacant Lot Locations", size(*1.2)) ///
subtitle ("Citywide RCT", size (*1.0))


But when I use geo2xy and project to mercator they no longer overlap on the x-axis. The y-axis looks okay.

use "C~phillyblockgscoord.dta", clear
geo2xy _Y _X, replace proj(mercator)
save "C~phillyblockscoordmerc.dta"

use "C~buffercluster.dta", clear
*reproject to coordinates
geo2xy latitude_y longitude_x, replace project(mercator)
save "C~cluster_locations_prj.dta", replace


use "C~phillyblockgsdb.dta", clear
spmap using "C~phillyblockscoordmerc.dta", id(id) ///
fcolor(eggshell) ocolor(grey) osize(thin) ///
point(data("C~cluster_locations_prj.dta") x(longitude_x) y(latitude_y) ///
size(*0.5) fcolor(black blue green) ocolor(white) ///
by (arm) ///
osize(vvthin)) ///
title ("Vacant Lot Locations", size(*1.2)) ///
subtitle ("Citywide RCT", size (*1.0))

[Robert Picard]

Most projections have parameters that determine the value of the projected values and in many case affect the look of the map (based on the projection's "point of view"). To keep things simple, geo2xy will calculate defaults based on the bounds of the coordinates to project.

You are using a mercator projection with default parameters. In that case, geo2xy uses the mid-longitude of the set of coordinates it is projecting. Here's an illustration of the issue where I overlay the location of Philadelphia and Washington, DC over a map of the lower 48 states:

clear
input _ID double(_Y _X)
1 39.951901 -75.166790
2 38.889696 -77.035258
end
save "dc_philly.dta", replace

geo2xy _Y _X, gen(latitude_y longitude_x) project(mercator)
ret list
save "dc_philly_proj.dta", replace

use "geo2xy_us_coor.dta", clear
drop if inlist(_ID, 14, 39, 42) // Alaska, Puerto Rico, Hawaii
geo2xy _Y _X, replace proj(mercator)
ret list
save "geo2xy_us_coor_proj.dta", replace

use "geo2xy_us_data.dta", clear
spmap using "geo2xy_us_coor_proj.dta", id(_ID) ///
fcolor(eggshell) ocolor(grey) osize(thin) ///
point(data("dc_philly_proj.dta") x(longitude_x) y(latitude_y))

graph export bad.png, width(800) replace

Here are the projection parameters used to project the cities:

. ret list

macros:
             r(aspect) : ".7329831224579442"
              r(model) : "Ellipsoid (6378137,298.257223563)"
              r(pname) : "Mercator"
               r(lon0) : "-76.101024"
                  r(f) : "298.257223563"
                  r(a) : "6378137"

and the one for the lower 48:

. ret list

macros:
             r(aspect) : ".5459560769389972"
              r(model) : "Ellipsoid (6378137,298.257223563)"
              r(pname) : "Mercator"
               r(lon0) : "-95.837867"
                  r(f) : "298.257223563"
                  r(a) : "6378137"

Since the lon0 parameter differs, the coordinates are not consistent between the two sets of coordinates and the map is wrong:

You could switch to the default projection (web_mercator) as it does not have any calculated parameters or you can use the parameters from the first projection on the second projection:

clear
input _ID double(_Y _X)
1 39.951901 -75.166790
2 38.889696 -77.035258
end
save "dc_philly.dta", replace

use "geo2xy_us_coor.dta", clear
drop if inlist(_ID, 14, 39, 42) // Alaska, Puerto Rico, Hawaii
geo2xy _Y _X, replace   proj(mercator)
ret list
local lon0 `r(lon0)'
local f `r(f)'
local a `r(a)'
save "geo2xy_us_coor_proj.dta", replace

use "dc_philly.dta", clear
geo2xy _Y _X, gen(latitude_y longitude_x)  proj(mercator, `a' `f' `lon0')
save "dc_philly_proj.dta", replace

use "geo2xy_us_data.dta", clear
spmap using "geo2xy_us_coor_proj.dta", id(_ID) ///
fcolor(eggshell) ocolor(grey) osize(thin) ///
point(data("dc_philly_proj.dta") x(longitude_x) y(latitude_y))

graph export mercator.png, width(800) replace

and here's the map generated by the above:

[John MacDonald]

Great that explains it! Thanks Robert for creating an excellent program.