#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
__author__ = "John Wieczorek"
__copyright__ = "Rauthiflor LLC"
__version__ = "datumshiftproj.py 2019-11-20T15:22-03:00"
__title__ = "Global estimates of worst-case datum shifts from WGS84"
# Use the Point class to determine haversine distances between coordinates
from point import *
# Use the PROJ4 python package for datum transformations
import pyproj
from pyproj import CRS
from pyproj import Transformer
from pyproj.enums import PJType
GRIDDEGREES=5
LNGMIN=-180
LNGMAX=180
LATMIN=-90
LATMAX=90.01
griddiffs=[]
gridpoints=[]
# print only the datum shifts for each cell, one per line
def printdiffsonly():
for diff in griddiffs:
print('%s' % (diff))
# print the datum shifts for each cell with the longitude and latitude of the SW corner of
# cell, one per line
def printvallnglat():
s = 'Value,Longitude,Latitude'
print(s)
i=0
for diff in griddiffs:
print('%s,%s,%s' % (diff, gridpoints[i].lng, gridpoints[i].lat))
i+=1
# True is string can be used to represent a floating point number
def is_float(string):
try:
float(string)
return True
except ValueError: # String is not a number
return False
# Similar to range(), but to allow a floating point number in a ranged for loop
def frange(start, stop=None, step=None):
# if stop and step argument is null set start=0.0 and step = 1.0
if stop == None:
stop = start + 0.0
start = 0.0
if step == None:
step = 1.0
while True:
if step > 0 and start >= stop:
break
elif step < 0 and start <= stop:
break
yield ("%g" % start) # return float number
start = start + step
def main():
'''
Create a csv file of datumshift maxima for ever grid cell of dimensions GRIDDEGREES
within:
LNGMIN, LNGMAX, LATMIN, LATMAX
Invoke without parameters as:
python datumshiftproj.py >> destfile.csv
'''
# Increment cells from LATMIN to LATMAX within LNGMIN to LNGMAX
for x in frange(LNGMIN,LNGMAX,GRIDDEGREES):
for y in frange(LATMIN,LATMAX,GRIDDEGREES):
# Create the cell with a diff of 0 and the coordinates of the SW corner of
# the cell
griddiffs.append(0)
gridpoints.append(Point(x,y))
# Prepare variables for maximum shift searches
# maxdiff - the overall biggest datum shift
# maxdiffpoint - the coordinates of the SW corner of the cell with the biggest
# datum shift
# maxcrs - the coordinate reference system of the SW corner of the cell with the
# biggest datum shift
maxdiff = 0
maxdiffpoint = None
maxcrs = None
# Cycle through all EPSG code for which the type is a geographic coordinate system
# and calculate datum shift for those that, like WGS84, use the Greenwich Meridian.
for code in pyproj.get_codes('EPSG',PJType.GEOGRAPHIC_CRS):
crs=CRS.from_user_input('epsg:%s' % code)
if crs.prime_meridian.name=='Greenwich':
i=0
t = Transformer.from_crs("epsg:%s" % code, 'epsg:4326', always_xy=True)
ops=t.operations
# For each grid cell in the current coordinate reference system, find the
# transformed corner coordinate in WGS84
for lng in frange(LNGMIN,LNGMAX,GRIDDEGREES):
for lat in frange(LATMIN,LATMAX,GRIDDEGREES):
p = Point(float(lng), float(lat))
x, y = t.transform(lng,lat)
p_wgs84 = Point(x,y)
try:
# Find the haversine distance between the point in the current
# coordinate reference system and the location with the same
# coordinate in WGS84
diff = p.haversine_distance(p_wgs84)
# If the datum shift for the current coordinate reference system
# is greater than any other for the same location so far, store it
if diff > griddiffs[i]:
griddiffs[i] = round(diff)
# If the datum shift for the current coordinate reference system
# is greater than any other so far, store it, along with the
# EPSG code for the coordinate reference system
if diff > maxdiff:
maxdiff=diff
maxdiffpoint = p
maxdiffwgs84 = Point(x,y)
maxcrs=crs
except ValueError:
pass
i+=1
# print(griddiffs)
# print("maxcrs: %s op:%s %.0fm at %s: WGS84: %s" % (maxcrs.to_wkt(), ops, maxdiff, maxdiffpoint, maxdiffwgs84))
# printvallnglat()
# print only the datum shifts for each cell
printdiffsonly()
if __name__ == "__main__":
main()