3

I've got a custom projection file that's applied to all data (per request of a client), and another team is requesting ids/coordinate pairs for assets. I tried ripping out X/Y off the point data, but I noticed the coordinates make no sense - they are not standard UTM northings/eastings.

For instance:

A point with real-world coords of 43.814794N, -76.02576907720452 has these coordinates in the weird custom projection: X: -534643.85116969, Y: 303444.080342568

The GCS is WGS84, but none of the online converters or formulas I can find work, and nobody has any clear instructions on how to create a formula. What I need is to write a Python function that will convert the eastings/northings to decimal degrees. Anyone have knowledge of how to develop such a formula?

Here's the info on the projection:

  False_Easting       |  0.000000000000000000

  False_Northing      |  0.000000000000000000

  Central_Meridian    |  -74.000000000000000000

  Standard_Parallel_1 | 42.000000000000000000

  Standard_Parallel_2 | 44.000000000000000000

  Latitude_Of_Origin  | 43.000000000000000000

And here's the info on the GCS (WGS84):

Name: GCS_WGS_1984

Angular Unit: Degree (0.0174532925199433)

Prime Meridian: Greenwich (0.0)

Datum: D_WGS_1984

  Spheroid: WGS_1984

    Semimajor Axis: 6378137.0

    Semiminor Axis: 6356752.314245179

    Inverse Flattening: 298.257223563
  • You tagged arcpy which works with ArcGIS and you have an Esri-style coordinate system definition. Do you have access to ArcGIS? – mkennedy Apr 17 at 18:03
  • 1
    Your coordinate reference system information is missing the projection algorithm. The parameters are for a conic projection--is it Lambert conformal conic or Albers equal area conic? – mkennedy Apr 17 at 18:04
  • @mkennedy yes, but I'm trying to avoid GP tools. I'd rather do calculations on X/Y extracted from Point objects on the fly in my loop that collects coordinate data. And neither - it is a custom conic projection. The client called it "Equidistant Conic." I posted everything available in the Proj file properties in ArcCatalog. – flintlockspecial Apr 17 at 18:06
  • 2
    Then you'll need to implement the equations on p114 of John P. Snyder's Map Projections: A Working Manual or try using pyproj maybe. – mkennedy Apr 17 at 18:13
3

If you decide to use pyproj (an excellent recommendation from @mkennedy), here is a projection string that may work:

+proj=eqdc +lat_0=43 +lat_1=42 +lat_2=44 +lon_0=-74 +datum=WGS84 +units=us-ft +no_defs

Using this online converter, it gives a result close to your example.

Input: +proj=longlat +datum=WGS84 +no_defs

Output: +proj=eqdc +lat_0=43 +lat_1=42 +lat_2=44 +lon_0=-74 +datum=WGS84 +units=us-ft +no_defs

43.814794, -76.02576907720452

Switch X <--> Y

Gives:

-534643.334276 ; 303440.613642

If you search this site, there are many excellent answers with examples of using pyproj; here is just one.

  • Thank you @Richard Morgan! I haven't had a chance to apply this yet (I work a different project Thurs-Fri), but I'll give it a shot Monday when I'm back at the client site! – flintlockspecial Apr 18 at 20:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.