Thank you Randy. It worked!
An alternative perhaps:
<SPAN class="keyword token">import</SPAN> arcpy x <SPAN class="operator token">=</SPAN> <SPAN class="number token">7079975.52661</SPAN> y <SPAN class="operator token">=</SPAN> <SPAN class="number token">2184310.50409</SPAN> <SPAN class="comment token"># WGS 1984 : (4326) Lat/Lon </SPAN> <SPAN class="comment token"># NAD_1983_2011_StatePlane_California_II_FIPS_0402_Ft_US - WKID: 6418 </SPAN> ptGeometry <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>PointGeometry<SPAN class="punctuation token">(</SPAN>arcpy<SPAN class="punctuation token">.</SPAN>Point<SPAN class="punctuation token">(</SPAN>x<SPAN class="punctuation token">,</SPAN>y<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">,</SPAN>arcpy<SPAN class="punctuation token">.</SPAN>SpatialReference<SPAN class="punctuation token">(</SPAN><SPAN class="number token">6418</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">.</SPAN>projectAs<SPAN class="punctuation token">(</SPAN>arcpy<SPAN class="punctuation token">.</SPAN>SpatialReference<SPAN class="punctuation token">(</SPAN><SPAN class="number token">4326</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN> <SPAN class="comment token"># print ptGeometry.JSON </SPAN> <SPAN class="keyword token">print</SPAN> ptGeometry<SPAN class="punctuation token">.</SPAN>firstPoint<SPAN class="punctuation token">.</SPAN>X<SPAN class="punctuation token">,</SPAN> ptGeometry<SPAN class="punctuation token">.</SPAN>firstPoint<SPAN class="punctuation token">.</SPAN>Y <SPAN class="comment token"># Result: -120.172262515 39.1457869508</SPAN><SPAN class="line-numbers-rows"><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN><SPAN></SPAN></SPAN>
Hi Chris,
I need to do similar conversion of California II State Plane in US Survey feet to geographic coordinates in decimal degrees. I am not totally comfortable with python but need to have it in python. I updated the parameters in California values. But not getting desired output. Here is the updated part of the script your wrote in California values:
A = 20925604.47 ## major radius of GRS 1980 ellipsoid, feet Ec = 0.0818191910435 ## eccentricity of GRD 1980 ellipsoid Ec2 = Ec * Ec ## eccentricity squared AngRad = 0.017453292519943299 ## number of radians in a degree Pi4 = 3.141592653582 / 4 ## Pi / 4 P1 = 38.33333333 * AngRad ## latitude of first standard parallel P2 = 39.83333333 * AngRad ## latitude of second standard parallel P0 = 37.66666667 * AngRad ## latitude of origin M0 = -122 * AngRad ## central meridian X0 = 6561666.66666667 ## False easting of central meridian, map units
And here is the snapshot of Data Source:
Unsure where the problem is. If it is due to no parameter for False northing or...?
Basically, the sample CA State Plane II x,y: 7079975.52661,2184310.50409 should give me: 39.1457869, -120.1722625 and if I run your code with above updated parameters, I am getting: 43.642324406, -120.048031426.
Any ideas how to fix that would be great!!!
def Ft_DD(x,y): ## This function converts coordinates from Florida State Plane North Zone NAD83 feet ## to Geographic coordinates in decimal degrees ## ## The formula this program is based on are from "Map Projections, ## A Working Manual" by John P. Snyder, U.S. Geological Survey ## Professional Paper 1395, 1987, pages 295-298 ## ##Function originally provided by Dan Saul, Washington State Department of Ecology ## ## Set up the coordinate system parameters. A = 20925604.47 ## major radius of GRS 1980 ellipsoid, feet Ec = 0.0818191910435 ## eccentricity of GRD 1980 ellipsoid Ec2 = Ec * Ec ## eccentricity squared AngRad = 0.017453292519943299 ## number of radians in a degree Pi4 = 3.141592653582 / 4 ## Pi / 4 P1 = 29.583333333333329000 * AngRad ## latitude of first standard parallel P2 = 30.750000000000000000 * AngRad ## latitude of second standard parallel P0 = 29.000000000000000000 * AngRad ## latitude of origin M0 = -84.5 * AngRad ## central meridian X0 = 1968500.0 ## False easting of central meridian, map units ## Calculate the coordinate system constants. m1 = math.cos(P1) / math.sqrt(1 - (Ec2 * pow((math.sin(P1)),2))) m2 = math.cos(P2) / math.sqrt(1 - (Ec2 * pow((math.sin(P2)),2))) t1 = math.tan(Pi4 - (P1 / 2)) / pow((1 - Ec * math.sin(P1)) / (1 + Ec * math.sin(P1)),(Ec/2)) t2 = math.tan(Pi4 - (P2 / 2)) / pow((1 - Ec * math.sin(P2)) / (1 + Ec * math.sin(P2)),(Ec/2)) t0 = math.tan(Pi4 - (P0 / 2)) / pow((1 - Ec * math.sin(P0)) / (1 + Ec * math.sin(P0)),(Ec/2)) n = math.log(m1 / m2) / math.log(t1 / t2) F = m1 / (n * pow(t1,n)) rho0 = A * F * pow(t0,n) ## Convert the coordinate to Latitude/Longitude. ## Calculate the Longitude. x = x - X0 Pi2 = Pi4 * 2 rho = math.sqrt(pow(x,2) + (pow(rho0 - y,2))) theta = math.atan(x/(rho0 - y)) t = pow(rho / (A * F),(1 / n)) LonR = theta / n + M0 x = x + X0 ## Estimate the Latitude Lat0 = Pi2 - (2 * math.atan(t)) ## Substitute the estimate into the iterative calculation that ## converges on the correct Latitude value. part1 = (1 - (Ec * math.sin(Lat0))) / (1 + (Ec * math.sin(Lat0))) LatR = Pi2 - (2 * math.atan(t * pow(part1,(Ec/2)))) Lat0 = LatR part1 = (1 - (Ec * math.sin(Lat0))) / (1 + (Ec * math.sin(Lat0))) ##LatR = Pi2 - (2 * math.atan2((t * (part1**(Ec/2))),1)) LatR = Pi2 - (2 * math.atan(t * pow(part1,(Ec/2)))) print LatR,LonR if (abs(LatR - Lat0) > 0.000000002): ## Convert from radians to degrees. Lat = LatR / AngRad Lon = LonR / AngRad print str(Lon) + "|" + str(Lat)
Angemeldete Mitglieder können Beiträge verfassen, Updates folgen und mehr. Neu hier? Registriere ein kostenloses Konto.
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.