diff options
Diffstat (limited to 'grid_tools/ntv2_to_gtiff.py')
| -rwxr-xr-x | grid_tools/ntv2_to_gtiff.py | 517 |
1 files changed, 517 insertions, 0 deletions
diff --git a/grid_tools/ntv2_to_gtiff.py b/grid_tools/ntv2_to_gtiff.py new file mode 100755 index 0000000..5c2449d --- /dev/null +++ b/grid_tools/ntv2_to_gtiff.py @@ -0,0 +1,517 @@ +#!/usr/bin/env python +############################################################################### +# $Id$ +# +# Project: PROJ +# Purpose: Convert a NTv2 file into an optimized GTiff +# Author: Even Rouault <even.rouault at spatialys.com> +# +############################################################################### +# Copyright (c) 2019, Even Rouault <even.rouault at spatialys.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. +############################################################################### + +# python ./ntv2_to_gtiff.py --copyright "Derived from work by IGN France. Open License https://www.etalab.gouv.fr/wp-content/uploads/2014/05/Open_Licence.pdf" --source-crs EPSG:4275 --target-crs EPSG:4171 /home/even/proj/proj-datumgrid/ntf_r93.gsb ntf_r93.tif +# python ./ntv2_to_gtiff.py --copyright "Derived from work by LGL-BW.DE. Data licence Germany - attribution - version 2.0: https://www.govdata.de/dl-de/by-2-0" --source-crs EPSG:4314 --target-crs EPSG:4258 /home/even/proj/proj-datumgrid/europe/BWTA2017.gsb BWTA2017.tif --do-not-write-accuracy-samples +# python ./ntv2_to_gtiff.py --copyright "Derived from work by Natural Resources Canada. Open Government Licence - Canada: http://open.canada.ca/en/open-government-licence-canada" --source-crs EPSG:4267 --target-crs EPSG:4269 /home/even/proj/proj-datumgrid/north-america/ntv2_0.gsb ntv2_0.tif +# python ./ntv2_to_gtiff.py --copyright "Derived from work by ICSM. Creative Commons Attribution 4.0: https://creativecommons.org/licenses/by/4.0/" --source-crs EPSG:4283 --target-crs EPSG:7844 /home/even/proj/proj-datumgrid/oceania/GDA94_GDA2020_conformal.gsb GDA94_GDA2020_conformal.tif + + +from osgeo import gdal +from osgeo import osr +from cloud_optimize_gtiff import generate_optimized_file +import argparse +import datetime +import math +import os +import struct + + +def get_args(): + parser = argparse.ArgumentParser( + description='Convert NTv2 grid into PROJ GeoTIFF.') + parser.add_argument('source', + help='Source NTv2 file') + parser.add_argument('dest', + help='Destination GeoTIFF file') + + parser.add_argument('--source-crs', dest='source_crs', required=True, + help='Source CRS as EPSG:XXXX or WKT') + + parser.add_argument('--target-crs', dest='target_crs', required=True, + help='Target CRS as EPSG:XXXX or WKT') + + parser.add_argument('--copyright', dest='copyright', required=True, + help='Copyright info') + + parser.add_argument('--description', dest='description', + help='Description') + + parser.add_argument('--do-not-write-accuracy-samples', dest='do_not_write_accuracy_samples', + action='store_true') + + parser.add_argument('--positive-longitude-shift-value', + dest='positive_longitude_shift_value', + choices=['east', 'west'], + default='east', + help='Whether positive values in the longitude_offset channel should be shift to the east or the west') + + parser.add_argument('--uint16-encoding', dest='uint16_encoding', + action='store_true', + help='Use uint16 storage with linear scaling/offseting') + + parser.add_argument('--datetime', dest='datetime', + help='Value for TIFF DateTime tag as YYYY:MM:DD HH:MM:SS, or "NONE" to not write it. If not specified, current date time is used') + + parser.add_argument('--accuracy-unit', dest='accuracy_unit', + choices=['arc-second', 'metre', 'unknown'], + help='Unit of accuracy channels') + + parser.add_argument('--area-of-use', dest='area_of_use', + help='Area of use') + + return parser.parse_args() + + +def get_year_month_day(src_date, src_basename): + if len(src_date) == 4 and src_basename == 'GS7783.GSB': + # CREATED=1991 + return int(src_date), 1, 1 + + if len(src_date) == 7 and src_basename == 'NB2783v2.GSB': + # CREATED=06/2011 + month = int(src_date[0:2]) + year = int(src_date[3:7]) + return year, month, 1 + + assert len(src_date) == 8 + if (src_date[2] == '-' and src_date[5] == '-') or \ + (src_date[2] == '/' and src_date[5] == '/'): + if src_basename.startswith('rdtrans') or \ + src_basename.startswith('rdcorr') or \ + src_basename.startswith('ntf_r93') or \ + src_basename.startswith('BWTA2017') or \ + src_basename.startswith('BETA2007') or \ + src_basename.startswith('D73_ETRS89_geo') or \ + src_basename.startswith('DLx_ETRS89_geo'): + # rdtrans2018.gsb has 22-11-18 & + # ntf_r93.gsb has 31/10/07, hence D-M-Y + day = int(src_date[0:2]) + month = int(src_date[3:5]) + year = int(src_date[6:8]) + else: + # CHENyx06a.gsb has 09-07-22 & ntv2_0.gsb and + # (other NRCan datasets) has 95-06-30, hence Y-M-D + year = int(src_date[0:2]) + month = int(src_date[3:5]) + day = int(src_date[6:8]) + if year >= 90: + year += 1900 + else: + assert year <= 50 + year += 2000 + else: + if src_basename in ('nzgd2kgrid0005.gsb', + 'A66_National_13_09_01.gsb', + 'National_84_02_07_01.gsb', + 'AT_GIS_GRID.gsb', + '100800401.gsb') or \ + src_basename.startswith('GDA94_GDA2020'): + # nzgd2kgrid0005 has 20111999, hence D-M-Y + day = int(src_date[0:2]) + month = int(src_date[2:4]) + year = int(src_date[4:8]) + elif src_basename == 'bd72lb72_etrs89lb08.gsb': + # bd72lb72_etrs89lb08 has 20142308, hence Y-D-M + year = int(src_date[0:4]) + day = int(src_date[4:6]) + month = int(src_date[6:8]) + else: + year = int(src_date[0:4]) + month = int(src_date[4:6]) + day = int(src_date[6:8]) + return year, month, day + + +def create_unoptimized_file(sourcefilename, tmpfilename, args): + src_ds = gdal.Open(sourcefilename) + assert src_ds.GetDriver().ShortName in ('NTv2', 'NTv1', 'CTable2') + subdatsets = [(sourcefilename, None)] + subdatsets += src_ds.GetSubDatasets() + + if src_ds.GetDriver().ShortName in ('NTv1', 'CTable2'): + args.do_not_write_accuracy_samples = True + + # Build a subgrids dict whose key is a parent grid name and the + # value the list of subgrids + subgrids = {} + for subds in subdatsets: + src_ds = gdal.Open(subds[0]) + parent_grid_name = src_ds.GetMetadataItem('PARENT') + if parent_grid_name is None or parent_grid_name == 'NONE': + continue + grid_name = src_ds.GetMetadataItem('SUB_NAME') + if parent_grid_name in subgrids: + subgrids[parent_grid_name].append(grid_name) + else: + subgrids[parent_grid_name] = [grid_name] + + src_basename = os.path.basename(args.source) + + if args.do_not_write_accuracy_samples: + nbands = 2 + else: + _, max = src_ds.GetRasterBand(3).ComputeRasterMinMax() + if max <= 0: + print('Omitting accuracy bands which contain only <= 0 values') + nbands = 2 + args.do_not_write_accuracy_samples = True + else: + if not args.accuracy_unit: + if src_basename in ('rdtrans2008.gsb', + 'rdtrans2018.gsb', + 'bd72lb72_etrs89lb08.gsb', + 'ntv2_0.gsb', + 'MAY76V20.gsb', + 'ABCSRSV4.GSB', + 'BC_27_05.GSB', + 'BC_93_05.GSB', + 'CQ77SCRS.GSB', + 'CRD27_00.GSB', + 'CRD93_00.GSB', + 'GS7783.GSB', + 'NA27SCRS.GSB', + 'NA83SCRS.GSB', + 'NB2783v2.GSB', + 'NB7783v2.GSB', + 'NS778302.GSB', + 'NVI93_05.GSB', + 'ON27CSv1.GSB', + 'ON76CSv1.GSB', + 'ON83CSv1.GSB', + 'PE7783V2.GSB', + 'SK27-98.GSB', + 'SK83-98.GSB', + 'TO27CSv1.GSB', + 'rdcorr2018.gsb',): + args.accuracy_unit = 'metre' + + elif src_basename in ('ntf_r93.gsb', + 'nzgd2kgrid0005.gsb', + 'OSTN15_NTv2_OSGBtoETRS.gsb', + 'CHENyx06a.gsb', + 'CHENyx06_ETRS.gsb', + 'A66_National_13_09_01.gsb', + 'National_84_02_07_01.gsb', + 'GDA94_GDA2020_conformal_and_distortion.gsb', + 'DLx_ETRS89_geo.gsb', + 'D73_ETRS89_geo.gsb'): + args.accuracy_unit = 'arc-second' + + else: + raise Exception( + '--accuracy-unit=arc-second/metre should be specified') + nbands = 4 + + compact_md = True if len(subdatsets) > 50 else False + + for idx_ifd, subds in enumerate(subdatsets): + src_ds = gdal.Open(subds[0]) + if src_ds.GetDriver().ShortName == 'NTv2': + assert src_ds.GetMetadataItem('GS_TYPE') == 'SECONDS' + + tmp_ds = gdal.GetDriverByName('GTiff').Create('/vsimem/tmp', + src_ds.RasterXSize, + src_ds.RasterYSize, + nbands, + gdal.GDT_Float32 if not args.uint16_encoding else gdal.GDT_UInt16) + src_crs = osr.SpatialReference() + src_crs.SetFromUserInput(args.source_crs) + tmp_ds.SetSpatialRef(src_crs) + tmp_ds.SetGeoTransform(src_ds.GetGeoTransform()) + tmp_ds.SetMetadataItem('AREA_OR_POINT', 'Point') + + grid_name = src_ds.GetMetadataItem('SUB_NAME') + if grid_name: + if src_basename == 'NVI93_05.GSB' and grid_name == 'NVIsib': + grid_name = grid_name + str(idx_ifd+1) + print('Fixing wrong SUB_NAME of NVI93_05.GSB to ' + grid_name) + tmp_ds.SetMetadataItem('grid_name', grid_name) + parent_grid_name = src_ds.GetMetadataItem('PARENT') + if parent_grid_name is None or parent_grid_name == 'NONE': + tmp_ds.SetMetadataItem('TYPE', 'HORIZONTAL_OFFSET') + else: + tmp_ds.SetMetadataItem('parent_grid_name', parent_grid_name) + if grid_name in subgrids: + tmp_ds.SetMetadataItem( + 'number_of_nested_grids', str(len(subgrids[grid_name]))) + + if idx_ifd == 0 or not compact_md: + # Indicates that positive shift values are corrections to the west ! + tmp_ds.GetRasterBand(2).SetMetadataItem('positive_value', + args.positive_longitude_shift_value) + + if args.uint16_encoding: + assert src_ds.GetDriver().ShortName == 'NTv2' + for i in (1, 2): + min, max = src_ds.GetRasterBand(i).ComputeRasterMinMax() + data = src_ds.GetRasterBand(i).ReadAsArray() + scale = (max - min) / 65535 + if i == 2 and args.positive_longitude_shift_value == 'east': + data = -data + data = (data - min) / scale + tmp_ds.GetRasterBand(i).WriteArray(data) + tmp_ds.GetRasterBand(i).SetOffset(min) + tmp_ds.GetRasterBand(i).SetScale(scale) + if idx_ifd == 0 or not compact_md: + tmp_ds.GetRasterBand(i).SetDescription( + 'latitude_offset' if i == 1 else 'longitude_offset') + tmp_ds.GetRasterBand(i).SetUnitType('arc-second') + + if nbands == 4: + for i in (3, 4): + min, max = src_ds.GetRasterBand(i).ComputeRasterMinMax() + data = src_ds.GetRasterBand(i).ReadAsArray() + scale = (max - min) / 65535 + if scale == 0: + data = 0 * data + else: + data = (data - min) / scale + tmp_ds.GetRasterBand(i).WriteArray(data) + tmp_ds.GetRasterBand(i).SetOffset(min) + tmp_ds.GetRasterBand(i).SetScale(scale) + if idx_ifd == 0 or not compact_md: + tmp_ds.GetRasterBand(i).SetDescription( + 'latitude_offset_accuracy' if i == 3 else 'longitude_offset_accuracy') + tmp_ds.GetRasterBand(i).SetUnitType(args.accuracy_unit) + + else: + for i in (1, 2): + data = src_ds.GetRasterBand(i).ReadRaster( + buf_type=gdal.GDT_Float32) + + if src_ds.GetDriver().ShortName == 'CTable2': + nvalues = src_ds.RasterXSize * src_ds.RasterYSize + out_data = b'' + # From radian to arc-seconds + for v in struct.unpack('f' * nvalues, data): + out_data += struct.pack('f', + v / math.pi * 180.0 * 3600) + data = out_data + + if i == 2 and args.positive_longitude_shift_value == 'east': + nvalues = src_ds.RasterXSize * src_ds.RasterYSize + out_data = b'' + for v in struct.unpack('f' * nvalues, data): + out_data += struct.pack('f', -v) + data = out_data + tmp_ds.GetRasterBand(i).WriteRaster(0, 0, src_ds.RasterXSize, src_ds.RasterYSize, + data) + if idx_ifd == 0 or not compact_md: + tmp_ds.GetRasterBand(i).SetDescription( + 'latitude_offset' if i == 1 else 'longitude_offset') + tmp_ds.GetRasterBand(i).SetUnitType('arc-second') + + if nbands == 4: + for i in (3, 4): + data = src_ds.GetRasterBand(i).ReadRaster() + tmp_ds.GetRasterBand(i).WriteRaster(0, 0, src_ds.RasterXSize, src_ds.RasterYSize, + data) + if idx_ifd == 0 or not compact_md: + tmp_ds.GetRasterBand(i).SetDescription( + 'latitude_offset_accuracy' if i == 3 else 'longitude_offset_accuracy') + tmp_ds.GetRasterBand(i).SetUnitType(args.accuracy_unit) + + dst_crs = osr.SpatialReference() + dst_crs.SetFromUserInput(args.target_crs) + dst_auth_name = dst_crs.GetAuthorityName(None) + dst_auth_code = dst_crs.GetAuthorityCode(None) + if idx_ifd == 0 or not compact_md: + if dst_auth_name == 'EPSG' and dst_auth_code: + tmp_ds.SetMetadataItem('target_crs_epsg_code', dst_auth_code) + else: + tmp_ds.SetMetadataItem( + 'target_crs_wkt', dst_crs.ExportToWkt(['FORMAT=WKT2_2018'])) + + if idx_ifd == 0: + desc = args.description + if not desc: + src_auth_name = src_crs.GetAuthorityName(None) + src_auth_code = src_crs.GetAuthorityCode(None) + + desc = src_crs.GetName() + if src_auth_name and src_auth_code: + desc += ' (' + src_auth_name + ':' + src_auth_code + ')' + desc += ' to ' + desc += dst_crs.GetName() + if dst_auth_name and dst_auth_code: + desc += ' (' + dst_auth_name + ':' + dst_auth_code + ')' + desc += '. Converted from ' + desc += src_basename + + extra_info = [] + version = src_ds.GetMetadataItem('VERSION') + if version: + version = version.strip() + if version not in ('NTv2.0',): + extra_info.append('version ' + version) + + src_date = src_ds.GetMetadataItem('UPDATED') + created_date = src_ds.GetMetadataItem('CREATED') + if src_date: + src_date = src_date.strip() + if created_date: + created_date = created_date.strip() + if not src_date: + src_date = created_date + + if created_date and src_date and len(src_date) < len(created_date): + # SK27-98.GSB + # CREATED=00-02-04 + # UPDATED=0-15-00 + # SK83-98.GSB + # CREATED=98-12-18 + # UPDATED=0-15-06 + src_date = created_date + + if src_date: + year, month, day = get_year_month_day( + src_date, src_basename) + + # Various sanity checks + assert day >= 1 and day <= 31 + assert month >= 1 and month <= 12 + assert year >= 1980 + assert year <= datetime.datetime.now().year + # assume agencies only work monday to friday... + # except in Belgium where they work on sundays + # and in NZ on saturdays + if src_basename not in ('nzgd2kgrid0005.gsb', + 'bd72lb72_etrs89lb08.gsb', + 'GS7783.GSB', + 'NB2783v2.GSB', + 'ON27CSv1.GSB', + 'ON76CSv1.GSB', + ): + assert datetime.datetime( + year, month, day).weekday() <= 4 + + # Sanity check that creation_date <= last_updated_date + if created_date: + year_created, month_created, day_created = get_year_month_day( + created_date, src_basename) + assert year_created * 10000 + month_created * 100 + \ + day_created <= year * 10000 + month * 100 + day + + extra_info.append( + 'last updated on %04d-%02d-%02d' % (year, month, day)) + + if extra_info: + desc += ' (' + ', '.join(extra_info) + ')' + + tmp_ds.SetMetadataItem('TIFFTAG_IMAGEDESCRIPTION', desc) + if args.copyright: + tmp_ds.SetMetadataItem('TIFFTAG_COPYRIGHT', args.copyright) + if args.datetime and args.datetime != 'NONE': + tmp_ds.SetMetadataItem('TIFFTAG_DATETIME', args.datetime) + if args.area_of_use: + tmp_ds.SetMetadataItem('area_of_use', args.area_of_use) + + options = ['PHOTOMETRIC=MINISBLACK', + 'COMPRESS=DEFLATE', + 'PREDICTOR=3' if not args.uint16_encoding else 'PREDICTOR=2', + 'INTERLEAVE=BAND', + 'GEOTIFF_VERSION=1.1'] + if tmp_ds.RasterXSize > 256 and tmp_ds.RasterYSize > 256: + options.append('TILED=YES') + else: + options.append('BLOCKYSIZE=' + str(tmp_ds.RasterYSize)) + if gdal.VSIStatL(tmpfilename) is not None: + options.append('APPEND_SUBDATASET=YES') + + assert gdal.GetDriverByName('GTiff').CreateCopy(tmpfilename, tmp_ds, + options=options) + + +def check(sourcefilename, destfilename, args): + src_ds = gdal.Open(sourcefilename) + assert src_ds.GetDriver().ShortName in ('NTv2', 'NTv1', 'CTable2') + src_subdatsets = [(sourcefilename, None)] + src_subdatsets += src_ds.GetSubDatasets() + + dst_ds = gdal.Open(destfilename) + dst_subdatsets = dst_ds.GetSubDatasets() + if not dst_subdatsets: + dst_subdatsets = [(destfilename, None)] + + assert len(src_subdatsets) == len(dst_subdatsets) + for src_subds, dst_subds in zip(src_subdatsets, dst_subdatsets): + src_ds = gdal.Open(src_subds[0]) + dst_ds = gdal.Open(dst_subds[0]) + if not args.uint16_encoding: + for i in range(min(src_ds.RasterCount, dst_ds.RasterCount)): + data = src_ds.GetRasterBand( + i+1).ReadRaster(buf_type=gdal.GDT_Float32) + + if src_ds.GetDriver().ShortName == 'CTable2': + nvalues = src_ds.RasterXSize * src_ds.RasterYSize + out_data = b'' + # From radian to arc-seconds + for v in struct.unpack('f' * nvalues, data): + out_data += struct.pack('f', + v / math.pi * 180.0 * 3600) + data = out_data + + if i+1 == 2 and args.positive_longitude_shift_value == 'east': + nvalues = src_ds.RasterXSize * src_ds.RasterYSize + out_data = b'' + for v in struct.unpack('f' * nvalues, data): + out_data += struct.pack('f', -v) + data = out_data + assert dst_ds.GetRasterBand(i+1).ReadRaster() == data + else: + import numpy as np + for i in range(min(src_ds.RasterCount, dst_ds.RasterCount)): + src_data = src_ds.GetRasterBand(i+1).ReadAsArray() + dst_data = dst_ds.GetRasterBand(i+1).ReadAsArray() + offset = dst_ds.GetRasterBand(i+1).GetOffset() + scale = dst_ds.GetRasterBand(i+1).GetScale() + dst_data = dst_data * scale + offset + if i+1 == 2 and args.positive_longitude_shift_value == 'east': + dst_data = -dst_data + max_error = np.max(abs(dst_data - src_data)) + assert max_error <= 1.01 * scale / 2, (max_error, scale / 2) + + +if __name__ == '__main__': + + args = get_args() + + tmpfilename = args.dest + '.tmp' + gdal.Unlink(tmpfilename) + + if not args.datetime and args.datetime != 'NONE': + args.datetime = datetime.date.today().strftime("%Y:%m:%d %H:%M:%S") + + create_unoptimized_file(args.source, tmpfilename, args) + generate_optimized_file(tmpfilename, args.dest) + check(args.source, args.dest, args) + + gdal.Unlink(tmpfilename) |