diff options
Diffstat (limited to 'grid_tools/convert_gr3df97a.py')
| -rwxr-xr-x | grid_tools/convert_gr3df97a.py | 146 |
1 files changed, 146 insertions, 0 deletions
diff --git a/grid_tools/convert_gr3df97a.py b/grid_tools/convert_gr3df97a.py new file mode 100755 index 0000000..f70927d --- /dev/null +++ b/grid_tools/convert_gr3df97a.py @@ -0,0 +1,146 @@ +#!/usr/bin/env python +############################################################################### +# $Id$ +# +# Project: PROJ +# Purpose: Convert IGN France geocentric gr3df97a.txt grid +# 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. +############################################################################### + +from osgeo import gdal +from osgeo import osr +import datetime +import struct + +# https://geodesie.ign.fr/contenu/fichiers/gr3df97a.txt + +with open('gr3df97a.txt', 'rt') as f: + count_val_rows = 0 + for line_no, l in enumerate(f.readlines()): + if l[-1] == '\n': + l = l[0:-1] + if l[-1] == '\r': + l = l[0:-1] + if line_no == 0: + assert l.startswith(' GR3D ') + elif line_no == 1: + assert l.startswith(' GR3D1 ') + tokens = [] + for t in l.split(' '): + if t: + tokens.append(t) + assert len(tokens) == 7 + minx = float(tokens[1]) + maxx = float(tokens[2]) + miny = float(tokens[3]) + maxy = float(tokens[4]) + resx = float(tokens[5]) + resy = float(tokens[6]) + assert minx < maxx + assert miny < maxy + assert resx > 0 + assert resy > 0 + cols = 1 + (maxx - minx) / resx + assert cols - int(cols + 0.5) < 1e-10 + cols = int(cols + 0.5) + rows = 1 + (maxy - miny) / resy + assert rows - int(rows + 0.5) < 1e-10 + rows = int(rows + 0.5) + ds = gdal.GetDriverByName('GTiff').Create( + '/vsimem/gr3df97a.tif', cols, rows, 3, gdal.GDT_Float32) + ds.SetMetadataItem( + 'TIFFTAG_COPYRIGHT', 'Derived from work by IGN France. Open License https://www.etalab.gouv.fr/wp-content/uploads/2014/05/Open_Licence.pdf') + ds.SetMetadataItem('TIFFTAG_IMAGEDESCRIPTION', + 'Geocentric translation from NTF (IGNF:NTF) to RGF93 (EPSG:4964). Converted from gr3df97a.txt') + datetime = datetime.date.today().strftime("%Y:%m:%d %H:%M:%S") + ds.SetMetadataItem('TIFFTAG_DATETIME', datetime) + ds.SetMetadataItem('AREA_OR_POINT', 'Point') + ds.SetMetadataItem('TYPE', 'GEOCENTRIC_TRANSLATION') + ds.SetMetadataItem('area_of_use', 'France') + ds.SetMetadataItem('source_crs_wkt', """GEODCRS["NTF cartesiennes", + DATUM["Nouvelle Triangulation Francaise", + ELLIPSOID["Clarke 1880 (IGN)",6378249.2,293.466021293627, + LENGTHUNIT["metre",1]]], + PRIMEM["Greenwich",0, + ANGLEUNIT["degree",0.0174532925199433]], + CS[Cartesian,3], + AXIS["(X)",geocentricX, + ORDER[1], + LENGTHUNIT["metre",1]], + AXIS["(Y)",geocentricY, + ORDER[2], + LENGTHUNIT["metre",1]], + AXIS["(Z)",geocentricZ, + ORDER[3], + LENGTHUNIT["metre",1]], + USAGE[ + SCOPE["unknown"], + AREA["FRANCE METROPOLITAINE (CORSE COMPRISE)"], + BBOX[41,-5.5,52,10]], + ID["IGNF","NTF"]]""") + ds.SetMetadataItem('target_crs_epsg_code', '4964') + ds.SetGeoTransform( + [minx - resx/2, resx, 0, maxy + resy/2, 0, -resy]) + ds.GetRasterBand(1).SetUnitType('metre') + ds.GetRasterBand(1).SetDescription('x_translation') + ds.GetRasterBand(2).SetUnitType('metre') + ds.GetRasterBand(2).SetDescription('y_translation') + ds.GetRasterBand(3).SetUnitType('metre') + ds.GetRasterBand(3).SetDescription('z_translation') + sr = osr.SpatialReference() + sr.ImportFromEPSG(4171) # RGF93 + ds.SetProjection(sr.ExportToWkt()) + elif line_no == 2: + assert l.startswith(' GR3D2 ') + elif line_no == 3: + assert l.startswith(' GR3D3 ') + else: + assert l.startswith('00002 ') + tokens = [] + for t in l.split(' '): + if t: + tokens.append(t) + assert len(tokens) == 8 + lon = float(tokens[1]) + lat = float(tokens[2]) + dx = float(tokens[3]) + dy = float(tokens[4]) + dz = float(tokens[5]) + iy = (line_no - 4) % rows + ix = (line_no - 4) // rows + assert abs(lon - (minx + ix * resx) + ) < 1e-10, (line_no, lon, minx + ix * resx) + assert abs(lat - (miny + iy * resy) + ) < 1e-10, (line_no, lat, miny + iy * resy) + ds.GetRasterBand(1).WriteRaster( + ix, rows - iy - 1, 1, 1, struct.pack('f', dx)) + ds.GetRasterBand(2).WriteRaster( + ix, rows - iy - 1, 1, 1, struct.pack('f', dy)) + ds.GetRasterBand(3).WriteRaster( + ix, rows - iy - 1, 1, 1, struct.pack('f', dz)) + count_val_rows += 1 + + assert count_val_rows == rows * cols + gdal.GetDriverByName('GTiff').CreateCopy('gr3df97a.tif', ds, options=[ + 'COMPRESS=DEFLATE', 'PREDICTOR=3', 'INTERLEAVE=BAND']) |