# How to test polar stereographic projection for a specific longitude intersection?

I'm writing a library that tests grids defined by map projections for intersections, and am running into an issue implementing the test for polar stereographic and rectangular projections - 'rectangular' meaning a grid where the sides of the projected grid coincide with latitude/longitude lines (e.g. mercator).

Without getting into the math too much, the formula to convert a lat/long coordinate to a north polar stereographic grid coordinate is:

E = (1 + sin(DxDyLatitude)) * EarthRadius
R = E * cos(Latitude) / (1 + sin(Latitude))
X = (sin(Longitude - Orientation) * -R) / IncrementX
Y = (cos(Longitude - Orientation) * R) / IncrementY
• IncrementX/Y** - the distance in meters between grid units of X/Y
• DxDyLatitude** - the latitude at which IncrementX/Y are valid
• Orientation** - the longitude that coincides with the +Y axis of the stereographic projection
• X/Y - the grid coordinate of the lat/long point

(** constant for a given grid)

To determine whether the top or bottom latitude of the rectangular grid intersects the stereographic grid, fix either X or Y at a grid boundary (left/right or bottom/top) and solve for the other. If the other is within left/right or bottom/top and the longitude at the same point is within the rectangular grid, the grids can be said to intersect. This is done by re-arranging the above equations:

Longitude = asin((X * IncrementX) / -R) + Orientation
Longitude = acos((Y * IncrementY) / R) + Orientation

However, determining whether the left or right longitude of the rectangular grid intersect seems not as easy due to the equation for R. When re-arranging the equations to solve for latitude instead of longitude, I am stuck at:

cos(Latitude) / (1 + sin(Latitude)) = X / ((sin(Longitude - Orientation) * -E) / IncrementX)
cos(Latitude) / (1 + sin(Latitude)) = Y / ((cos(Longitude - Orientation) * E) / IncrementY)

I don't know of any trigonometric identity that reduces cos(x) / (1 + sin(x)), or if the equation can be re-arranged so that there is such an identity, or of any other mathematical equation for stereographic projections that will allow me to do this.

• I don't have an answer, but if you just want to do it use the transformation functions in proj4 or rgdal. – mdsumner Nov 10 '15 at 0:43