Monday, March 27, 2017

Convert WGS_1984 MSL to absolute altitude (AGL)

Here's the latest open source code coming out of the effort I've been supporting... the code converts WGS_1984 MSL to absolute altitude (AGL).  You can find the code on GitHub under the eleveation2agl project.
Determining AGL (elevation2agl.py)...

Example Inputs...
 WGS_1984 location                              : -74.5659708,39.4505594
 Provided reference elevation (WGS_1984 MSL)    : 1660 feet above MSL

Ground level calculations...
  NAD_1983 location (used w/ ground level query): -74.56596924794374,39.45055141406901
  Gathered ground level altitude                : 59.73 feet above MSL

AGL Calculations...
  Absolute altitude - above ground level (AGL)  : 1600.27 feet 

Here's a quicklook at the latest code/script at the time of this writing...

#!/usr/local/bin/python3
#
# Purpose: Converts WGS_1984 MSL to absolute altitude (AGL)
#
# Example usage:
#    python elevation2agl.py
#
# References:
# view-source:http://tagis.dep.wv.gov/convert/
# http://tagis.dep.wv.gov/arcgis/rest/services/Utilities/Geometry/GeometryServer/project
# https://developers.arcgis.com/javascript/3/jsapi/geometryservice-amd.html
# Projected Coordinate Systems https://developers.arcgis.com/javascript/3/jshelp/pcs.html
#
# https://tasks.arcgisonline.com/arcgis/rest/services/Geometry/GeometryServer/
# https://tasks.arcgisonline.com/arcgis/rest/services/Geometry/GeometryServer/project?inSR=4326&outSR=104223&geometries=%7B%22geometryType%22%3A%22esriGeometryPoint%22%2C%22geometries%22%3A%5B%7B%22x%22%3A-117%2C%22y%22%3A34%7D%5D%7D&transformation=&transformForward=true&f=json
# Input Spatial Reference:    4326 (GCS_WGS_1984)
# Output Spatial Reference: 104223 (GCS_NAD_1983_CORS96)
# Geometries: {"geometryType":"esriGeometryPoint","geometries":[{"x":-117,"y":34}]}
#
# "Without digging into the details, the datum realizations (such as 104223 vs. 4269) are not mathematically identical,
# but for almost any real-world purpose they are so close it doesn't matter. What does matter is that the grid point
# spacing of the 1/3 arc-second data is about 10 meters, and this relatively low resolution elevation interpolation
# -- and other uncertainties inherent in the capture of these data -- overcomes the differences in the datum math models
# by at least two to three orders of magnitude.  See FAQ https://www2.usgs.gov/faq/categories/9865/3624 for more
# discussion on elevation accuracy of the EPQS"..."If you are attempting to improve the positional accuracy of the data,
# then yes, this conversion is unnecessary. Changes in positions will be small relative to the point spacing, and change
# in actual overall accuracy will be insignificant. But it doesn't do any harm, either, and sometimes conversions like
# this are desirable for other reasons...achieving consistency in metadata, or making data friendlier to particular
# software processes, for example." -- USGS TNM Service Desk
#
# ########################################################################################################################
__author__  = ["Robert James Liguori (HAI/Gliesian", "Lana Manovych (ATAC)"]
__script__  = "elevation2Agl.py"
__version__ = "0.1.10 03-27-17"
__since__   = "0.1.0  02-10-16"
__status__  = "Deployment"

import requests
import json

#######################
# USER INPUT (WGS 1984)
######################
referenceElevation = 1660
msg84_lon = str(-74.5659708)
msg84_lat = str(39.4505594)

#######################################################################################
# Convert from WGS84 to North American Datum 1983 (NAD83) using arcgis geometry service
# (REF: https://nationalmap.gov/epqs/)
#######################################################################################
INPUT_DATUM  = "GCS_WGS_1984:4326"
TARGET_DATUM = "GCS_NAD_1983_CORS96:104223"
payload = {'inSR': INPUT_DATUM.split(":")[1],
                          'outSR': TARGET_DATUM.split(":")[1],
                          'geometries': json.dumps({'geometryType': 'esriGeometryPoint', 'geometries': [{'x':msg84_lon,'y': msg84_lat }]}),
                          'transformation': '',
                          'transformForward': 'true',
                          'f': 'json'}
r = requests.put('https://tasks.arcgisonline.com/arcgis/rest/services/Geometry/GeometryServer/project',
                  params=payload)
data = json.loads (r.text)

for rows in data["geometries"]:
    nad83_lon = str(rows['x'])
    nad83_lat = str(rows['y'])

###############
# DETERMINE AGL
###############

# Ground level altitude/elevation above mean sea level
# e.g. http://nationalmap.gov/epqs/pqs.php?x=-74.5659708&y=39.4505594&units=Feet&output=xml
line = requests.get('http://nationalmap.gov/epqs/pqs.php?x=' + nad83_lon + '&y=' + nad83_lat + '&units=Feet&output=xml')
if line == -1000000:
    exit ("service unable to find data at the requested point")
groundLevelAltitude = (str(line.content).split('Elevation>'))[1].split('<')[0]

# Rotorcraft's absolute altitude above ground level (AGL)
absoluteAltitude = referenceElevation - float(groundLevelAltitude)

###################
# PRINT CALCULATION
###################
print("\nDetermining AGL (elevation2agl.py)...")

print("\nExample Inputs...")
print(" WGS_1984 location                              : " + msg84_lon +"," + msg84_lat)
print(" Provided reference elevation (WGS_1984 MSL)    : " + str(referenceElevation) + " feet above MSL")

print("\nGround level calculations...")
print("  NAD_1983 location (used w/ ground level query): " + nad83_lon +"," + nad83_lat)
print("  Gathered ground level altitude                : " + str(groundLevelAltitude) + " feet above MSL")

print("\nAGL Calculations...")
print("  Absolute altitude - above ground level (AGL)  : " + str(absoluteAltitude) + " feet above ground")

No comments:

Post a Comment

Landing Spicoli in the Wind

Many know by now that I have three Blade 200 S RC helicopters in my fleet: Speck - for training Spot - for maintenance exercises Spicol...