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
nmtoken
  • 13,355
  • 5
  • 38
  • 87
flintlockspecial
  • 306
  • 1
  • 8
  • 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 '19 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 '19 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 '19 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 '19 at 18:13

1 Answers1

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.

Community
  • 1
Richard Morgan
  • 2,819
  • 1
  • 21
  • 28
  • 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 '19 at 20:03