JezK
Edit File: crs2crs2grid.py
#!/usr/bin/env python3 # ****************************************************************************** # # Project: GDAL # Purpose: Use HTDP to generate PROJ.4 compatible datum grid shift files. # Author: Frank Warmerdam, warmerdam@pobox.com # # See also: http://www.ngs.noaa.gov/TOOLS/Htdp/Htdp.shtml # http://trac.osgeo.org/proj/wiki/HTDPGrids # # ****************************************************************************** # Copyright (c) 2012, Frank Warmerdam <warmerdam@pobox.com> # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included # in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER # DEALINGS IN THE SOFTWARE. # ****************************************************************************** import os import sys import numpy from osgeo import gdal, gdal_array # Input looks like this: """ PNT_0_0 LATITUDE 40 00 0.00000 N 40 00 0.02344 N 2.06 mm/yr north LONGITUDE 117 00 0.00000 W 117 00 0.03765 W -1.22 mm/yr east ELLIP. HT. 0.000 -0.607 m -1.34 mm/yr up X -2221242.768 -2221243.142 m -0.02 mm/yr Y -4359434.393 -4359433.159 m 2.64 mm/yr Z 4077985.572 4077985.735 m 0.72 mm/yr """ def next_point(fd): line = fd.readline().strip() while line != "" and line.find("PNT_") == -1: line = fd.readline() if line == "": return None name_tokens = line.split("_") # Read LATITUDE line line = fd.readline().strip() tokens = line.split() lat_src = float(tokens[1]) + float(tokens[2]) / 60.0 + float(tokens[3]) / 3600.0 lat_dst = float(tokens[5]) + float(tokens[6]) / 60.0 + float(tokens[7]) / 3600.0 line = fd.readline().strip() tokens = line.split() lon_src = float(tokens[1]) + float(tokens[2]) / 60.0 + float(tokens[3]) / 3600.0 lon_dst = float(tokens[5]) + float(tokens[6]) / 60.0 + float(tokens[7]) / 3600.0 return ( int(name_tokens[1]), int(name_tokens[2]), lat_src, lon_src, lat_dst, lon_dst, ) def read_grid_crs_to_crs(filename, shape): fd = open(filename) grid = numpy.zeros(shape) # report the file header defining the transformation for _ in range(5): print(fd.readline().rstrip()) points_found = 0 ptuple = next_point(fd) while ptuple is not None: grid[0, ptuple[1], ptuple[0]] = ptuple[4] - ptuple[2] grid[1, ptuple[1], ptuple[0]] = ptuple[5] - ptuple[3] points_found = points_found + 1 ptuple = next_point(fd) if points_found < shape[0] * shape[1]: print("points found: ", points_found) print("points expected:", shape[1] * shape[2]) return None return grid ############################################################################### # This function creates a regular grid of lat/long values with one # "band" for latitude, and one for longitude. # def new_create_grid(griddef): lon_start = -1 * griddef[0] lon_end = -1 * griddef[2] lat_start = griddef[1] lat_end = griddef[3] lon_steps = griddef[4] lat_steps = griddef[5] lat_axis = numpy.linspace(lat_start, lat_end, lat_steps) lon_axis = numpy.linspace(lon_start, lon_end, lon_steps) lon_list = [lon_axis] * lat_steps lon_band = numpy.vstack(lon_list) lat_list = [lat_axis] * lon_steps lat_band = numpy.column_stack(lat_list) return numpy.array([lon_band, lat_band]) ############################################################################## # This function writes a grid out in form suitable to use as input to the # htdp program. def write_grid(grid, out_filename): fd_out = open(out_filename, "w") for i in range(grid.shape[2]): for j in range(grid.shape[1]): fd_out.write('%f %f 0 "PNT_%d_%d"\n' % (grid[1, j, i], grid[0, j, i], i, j)) fd_out.close() ############################################################################## # Write the resulting grid out in GeoTIFF format. def write_gdal_grid(filename, grid, griddef): ps_x = (griddef[2] - griddef[0]) / (griddef[4] - 1) ps_y = (griddef[3] - griddef[1]) / (griddef[5] - 1) geotransform = ( griddef[0] - ps_x * 0.5, ps_x, 0.0, griddef[1] - ps_y * 0.5, 0.0, ps_y, ) grid = grid.astype(numpy.float32) ds = gdal_array.SaveArray(grid, filename, format="CTable2") ds.SetGeoTransform(geotransform) ############################################################################# def write_control( control_fn, out_grid_fn, in_grid_fn, src_crs_id, src_crs_date, dst_crs_id, dst_crs_date, ): # start_date, end_date should be something like "2011.0" control_template = """ 4 %s %d %d 2 %s 2 %s 3 %s 0 0 """ control_filled = control_template % ( out_grid_fn, src_crs_id, dst_crs_id, src_crs_date, dst_crs_date, in_grid_fn, ) open(control_fn, "w").write(control_filled) ############################################################################# def Usage(brief=1): print( """ Usage: crs2crs2grid.py <src_crs_id> <src_crs_date> <dst_crs_id> <dst_crs_year> [-griddef <ul_lon> <ul_lat> <ll_lon> <ll_lat> <lon_count> <lat_count>] [-htdp <path_to_exe>] [-wrkdir <dirpath>] [-kwf] -o <output_grid_name> -griddef: by default the following values for roughly the continental USA at a six minute step size are used: -127 50 -66 25 251 611 -kwf: keep working files in the working directory for review. eg. crs2crs2grid.py 29 2002.0 8 2002.0 -o nad83_2002.ct2 """ ) if brief == 0: print( """ The output file will be in CTable2 format suitable for use with PROJ.4 +nadgrids= directive. Format dates like 2002.0 (for the start of 2002) CRS Ids ------- 1...NAD_83(2011) (North America tectonic plate fixed) 29...NAD_83(CORS96) (NAD_83(2011) will be used) 30...NAD_83(2007) (NAD_83(2011) will be used) 2...NAD_83(PA11) (Pacific tectonic plate fixed) 31...NAD_83(PACP00) (NAD 83(PA11) will be used) 3...NAD_83(MA11) (Mariana tectonic plate fixed) 32...NAD_83(MARP00) (NAD_83(MA11) will be used) 4...WGS_72 16...ITRF92 5...WGS_84(transit) = NAD_83(2011) 17...ITRF93 6...WGS_84(G730) = ITRF92 18...ITRF94 = ITRF96 7...WGS_84(G873) = ITRF96 19...ITRF96 8...WGS_84(G1150) = ITRF2000 20...ITRF97 9...PNEOS_90 = ITRF90 21...IGS97 = ITRF97 10...NEOS_90 = ITRF90 22...ITRF2000 11...SIO/MIT_92 = ITRF91 23...IGS00 = ITRF2000 12...ITRF88 24...IGb00 = ITRF2000 13...ITRF89 25...ITRF2005 14...ITRF90 26...IGS05 = ITRF2005 15...ITRF91 27...ITRF2008 28...IGS08 = ITRF2008 """ ) return 2 ############################################################################# # Main def main(argv=sys.argv): # Default GDAL argument parsing. argv = gdal.GeneralCmdLineProcessor(argv) if argv is None: return 0 if len(argv) == 1: return Usage(brief=0) # Script argument defaults src_crs_id = None src_crs_date = None dst_crs_id = None dst_crs_date = None # Decent representation of continental US griddef = (-127.0, 50.0, -66.0, 25.0, 611, 251) htdp_path = "htdp" wrkdir = "." kwf = 0 output_grid_name = None # Script argument parsing. i = 1 while i < len(argv): if argv[i] == "-griddef" and i < len(argv) - 6: griddef = ( float(argv[i + 1]), float(argv[i + 2]), float(argv[i + 3]), float(argv[i + 4]), float(argv[i + 5]), float(argv[i + 6]), ) i = i + 6 elif argv[i] == "-htdp" and i < len(argv) - 1: htdp_path = argv[i + 1] i = i + 1 elif argv[i] == "-kwf": kwf = 1 elif argv[i] == "-wrkdir" and i < len(argv) - 1: wrkdir = argv[i + 1] i = i + 1 elif argv[i] == "-o" and i < len(argv) - 1: output_grid_name = argv[i + 1] i = i + 1 elif argv[i] == "-h" or argv[i] == "--help": return Usage(brief=0) elif argv[i][0] == "-": print("Urecognised argument: " + argv[i]) return Usage() elif src_crs_id is None: src_crs_id = int(argv[i]) elif src_crs_date is None: src_crs_date = argv[i] elif dst_crs_id is None: dst_crs_id = int(argv[i]) elif dst_crs_date is None: dst_crs_date = argv[i] else: print("Urecognised argument: " + argv[i]) return Usage() i = i + 1 # next argument if output_grid_name is None: print("Missing output grid name (-o)") return Usage() if dst_crs_date is None: print( "Source and Destination CRS Ids and Dates are mandatory, " "not all provided." ) return Usage() # Do a bit of validation of parameters. if src_crs_id < 1 or src_crs_id > 32 or dst_crs_id < 1 or dst_crs_id > 32: print( "Invalid source or destination CRS Id %d and %d." % (src_crs_id, dst_crs_id) ) return Usage(brief=0) if ( float(src_crs_date) < 1700.0 or float(src_crs_date) > 2300.0 or float(dst_crs_date) < 1700.0 or float(dst_crs_date) > 2300.0 ): print( "Source or destination CRS date seems odd %s and %s." % (src_crs_date, dst_crs_date) ) return Usage(brief=0) # Prepare out set of working file names. in_grid_fn = wrkdir + "/crs2crs_input.txt" out_grid_fn = wrkdir + "/crs2crs_output.txt" control_fn = wrkdir + "/crs2crs_control.txt" # Write out the source grid file. grid = new_create_grid(griddef) write_grid(grid, in_grid_fn) write_control( control_fn, out_grid_fn, in_grid_fn, src_crs_id, src_crs_date, dst_crs_id, dst_crs_date, ) # Run htdp to transform the data. try: os.unlink(out_grid_fn) except OSError: pass rc = os.system(htdp_path + " < " + control_fn) if rc != 0: print("htdp run failed!") return 1 print("htdp run complete.") adjustment = read_grid_crs_to_crs(out_grid_fn, grid.shape) if adjustment is None: return 1 # Convert shifts to radians adjustment = adjustment * (3.14159265358979323846 / 180.0) # write to output output grid file. write_gdal_grid(output_grid_name, adjustment, griddef) # cleanup working files unless they have been requested to remain. if kwf == 0: os.unlink(in_grid_fn) os.unlink(out_grid_fn) os.unlink(control_fn) print("Processing complete: see " + output_grid_name) return 0 if __name__ == "__main__": sys.exit(main(sys.argv))