From 57f9efad9a7e33722eed0b22ffe23c1df8a9652d Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 10 Mar 2020 21:34:08 +0100 Subject: Add nc_dittt/ grids for New Caledonia (fixes #14) --- CMakeLists.txt | 1 + agency.json | 6 ++ copyright_and_licenses.csv | 4 ++ grid_tools/convert_gr3dnc01b.py | 121 ----------------------------------- grid_tools/convert_gr3dnc03a.py | 121 ----------------------------------- nc_dittt/convert_gr3dnc01b.py | 138 ++++++++++++++++++++++++++++++++++++++++ nc_dittt/convert_gr3dnc02b.py | 138 ++++++++++++++++++++++++++++++++++++++++ nc_dittt/convert_gr3dnc03a.py | 138 ++++++++++++++++++++++++++++++++++++++++ nc_dittt/nc_dittt_README.txt | 33 ++++++++++ nc_dittt/nc_dittt_gr3dnc01b.tif | Bin 0 -> 62390 bytes nc_dittt/nc_dittt_gr3dnc02b.tif | Bin 0 -> 7975 bytes nc_dittt/nc_dittt_gr3dnc03a.tif | Bin 0 -> 3274 bytes travis/expected_main.lst | 4 ++ 13 files changed, 462 insertions(+), 242 deletions(-) delete mode 100755 grid_tools/convert_gr3dnc01b.py delete mode 100755 grid_tools/convert_gr3dnc03a.py create mode 100755 nc_dittt/convert_gr3dnc01b.py create mode 100755 nc_dittt/convert_gr3dnc02b.py create mode 100755 nc_dittt/convert_gr3dnc03a.py create mode 100644 nc_dittt/nc_dittt_README.txt create mode 100644 nc_dittt/nc_dittt_gr3dnc01b.tif create mode 100644 nc_dittt/nc_dittt_gr3dnc02b.tif create mode 100644 nc_dittt/nc_dittt_gr3dnc03a.tif diff --git a/CMakeLists.txt b/CMakeLists.txt index 7bc5983..289a62b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -33,6 +33,7 @@ set(CPACK_SOURCE_INSTALLED_DIRECTORIES "${CMAKE_CURRENT_SOURCE_DIR}/eur_nkg/" "." "${CMAKE_CURRENT_SOURCE_DIR}/fr_ign/" "." "${CMAKE_CURRENT_SOURCE_DIR}/is_lmi/" "." + "${CMAKE_CURRENT_SOURCE_DIR}/nc_dittt/" "." "${CMAKE_CURRENT_SOURCE_DIR}/nl_nsgi/" "." "${CMAKE_CURRENT_SOURCE_DIR}/nz_linz/" "." "${CMAKE_CURRENT_SOURCE_DIR}/pt_dgt/" "." diff --git a/agency.json b/agency.json index 80bf6df..a9cae5e 100644 --- a/agency.json +++ b/agency.json @@ -95,6 +95,12 @@ "agency":"National Land Survey of Iceland", "url": "https://atlas.lmi.is/LmiData/index.php?id=626468364600" }, + { + "id": "nc_dittt", + "country": "New Caledonia", + "agency":"Gouvernement de Nouvelle Calédonie - DITTT", + "url": "https://dittt.gouv.nc/circe-pour-la-nouvelle-caledonie" + }, { "id": "nl_nsgi", "country": "Netherlands", diff --git a/copyright_and_licenses.csv b/copyright_and_licenses.csv index 38cba73..6c5106d 100644 --- a/copyright_and_licenses.csv +++ b/copyright_and_licenses.csv @@ -137,6 +137,10 @@ is_lmi_ISN93_ISN2016.tif,"Copyright 2017-2019, National Land Survey of Iceland", is_lmi_ISN_vel_beta.tif,"Copyright 2017-2019, National Land Survey of Iceland",CC-BY-4.0 is_lmi_README.txt,Disclaimed,Public domain NKG,Nordic Geodetic Commission,CC-BY-4.0 +nc_dittt_gr3dnc01b.tif,Gouvernement de Nouvelle Calédonie - DITTT,Open License France - https://www.etalab.gouv.fr/wp-content/uploads/2014/05/Open_Licence.pdf +nc_dittt_gr3dnc02b.tif,Gouvernement de Nouvelle Calédonie - DITTT,Open License France - https://www.etalab.gouv.fr/wp-content/uploads/2014/05/Open_Licence.pdf +nc_dittt_gr3dnc03a.tif,Gouvernement de Nouvelle Calédonie - DITTT,Open License France - https://www.etalab.gouv.fr/wp-content/uploads/2014/05/Open_Licence.pdf +nc_dittt_README.txt,Disclaimed,Public domain nl_nsgi_naptrans2018.tif,Nederlandse Samenwerking Geodetische Infrastructuur (NSGI),CC-BY-4.0 nl_nsgi_nlgeo2018.tif,Nederlandse Samenwerking Geodetische Infrastructuur (NSGI),CC-BY-4.0 nl_nsgi_rdcorr2018.tif,Nederlandse Samenwerking Geodetische Infrastructuur (NSGI),CC-BY-4.0 diff --git a/grid_tools/convert_gr3dnc01b.py b/grid_tools/convert_gr3dnc01b.py deleted file mode 100755 index b3ed17e..0000000 --- a/grid_tools/convert_gr3dnc01b.py +++ /dev/null @@ -1,121 +0,0 @@ -#!/usr/bin/env python -############################################################################### -# $Id$ -# -# Project: PROJ -# Purpose: Convert New Caledonia geocentric gr3dnc01b.mnt grid -# Author: Even Rouault -# -############################################################################### -# Copyright (c) 2019, Even Rouault -# -# 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 - -with open('gr3dnc01b.mnt', 'rt') as f: - count_val_rows = 0 - line_no = 0 - for l in f.readlines(): - if l[-1] == '\n': - l = l[0:-1] - if l[-1] == '\r': - l = l[0:-1] - if not l: - continue - if line_no == 0: - tokens = [] - for t in l.split(' '): - if t: - tokens.append(t) - assert len(tokens) >= 13 - minx = float(tokens[0]) - maxx = float(tokens[1]) - miny = float(tokens[2]) - maxy = float(tokens[3]) - resx = float(tokens[4]) - resy = float(tokens[5]) - 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/gr3dnc01b.tif', cols, rows, 3, gdal.GDT_Float32) - ds.SetMetadataItem( - 'TIFFTAG_COPYRIGHT', 'Derived from work by Service Topographique, DITTT, GNC. License unclear') - ds.SetMetadataItem('TIFFTAG_IMAGEDESCRIPTION', - 'Geocentric translation from IGN72 GRANDE TERRE (EPSG:4662) to RGNC91-93 (EPSG:4749). Converted from gr3dnc01b.mnt') - 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_code', '4662') - ds.SetMetadataItem('target_crs_epsg_code', '4749') - 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(4749) - ds.SetProjection(sr.ExportToWkt()) - else: - tokens = [] - for t in l.split(' '): - if t: - tokens.append(t) - assert len(tokens) == 6 - lon = float(tokens[0]) - lat = float(tokens[1]) - dx = float(tokens[2]) - dy = float(tokens[3]) - dz = float(tokens[4]) - iy = (line_no - 1) % rows - ix = (line_no - 1) // 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 - - line_no += 1 - - assert count_val_rows == rows * cols - gdal.GetDriverByName('GTiff').CreateCopy('gr3dnc01b.tif', ds, options=[ - 'COMPRESS=DEFLATE', 'PREDICTOR=3', 'INTERLEAVE=BAND']) diff --git a/grid_tools/convert_gr3dnc03a.py b/grid_tools/convert_gr3dnc03a.py deleted file mode 100755 index 85c35a3..0000000 --- a/grid_tools/convert_gr3dnc03a.py +++ /dev/null @@ -1,121 +0,0 @@ -#!/usr/bin/env python -############################################################################### -# $Id$ -# -# Project: PROJ -# Purpose: Convert New Caledonia geocentric gr3dnc03a.mnt grid -# Author: Even Rouault -# -############################################################################### -# Copyright (c) 2019, Even Rouault -# -# 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 - -with open('gr3dnc03a.mnt', 'rt') as f: - count_val_rows = 0 - line_no = 0 - for l in f.readlines(): - if l[-1] == '\n': - l = l[0:-1] - if l[-1] == '\r': - l = l[0:-1] - if not l: - continue - if line_no == 0: - tokens = [] - for t in l.split(' '): - if t: - tokens.append(t) - assert len(tokens) >= 13 - minx = float(tokens[0]) - maxx = float(tokens[1]) - miny = float(tokens[2]) - maxy = float(tokens[3]) - resx = float(tokens[4]) - resy = float(tokens[5]) - 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/gr3dnc03a.tif', cols, rows, 3, gdal.GDT_Float32) - ds.SetMetadataItem( - 'TIFFTAG_COPYRIGHT', 'Derived from work by Service Topographique, DITTT, GNC. License unclear') - ds.SetMetadataItem('TIFFTAG_IMAGEDESCRIPTION', - 'Geocentric translation from NEA74 Noumea (EPSG:4644) to RGNC91-93 (EPSG:4749). Converted from gr3dnc03a.mnt') - 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_code', '4644') - ds.SetMetadataItem('target_crs_epsg_code', '4749') - 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(4749) - ds.SetProjection(sr.ExportToWkt()) - else: - tokens = [] - for t in l.split(' '): - if t: - tokens.append(t) - assert len(tokens) == 6 - lon = float(tokens[0]) - lat = float(tokens[1]) - dx = float(tokens[2]) - dy = float(tokens[3]) - dz = float(tokens[4]) - iy = (line_no - 1) % rows - ix = (line_no - 1) // 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 - - line_no += 1 - - assert count_val_rows == rows * cols - gdal.GetDriverByName('GTiff').CreateCopy('gr3dnc03a.tif', ds, options=[ - 'COMPRESS=DEFLATE', 'PREDICTOR=3', 'INTERLEAVE=BAND']) diff --git a/nc_dittt/convert_gr3dnc01b.py b/nc_dittt/convert_gr3dnc01b.py new file mode 100755 index 0000000..c169de0 --- /dev/null +++ b/nc_dittt/convert_gr3dnc01b.py @@ -0,0 +1,138 @@ +#!/usr/bin/env python +############################################################################### +# $Id$ +# +# Project: PROJ +# Purpose: Convert New Caledonia geocentric gr3dnc01b.mnt grid +# Author: Even Rouault +# +############################################################################### +# Copyright (c) 2019, Even Rouault +# +# 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 + +with open('gr3dnc01b.mnt', 'rt') as f: + count_val_rows = 0 + line_no = 0 + for l in f.readlines(): + if l[-1] == '\n': + l = l[0:-1] + if l[-1] == '\r': + l = l[0:-1] + if not l: + continue + if line_no == 0: + tokens = [] + for t in l.split(' '): + if t: + tokens.append(t) + assert len(tokens) >= 13 + minx = float(tokens[0]) + maxx = float(tokens[1]) + miny = float(tokens[2]) + maxy = float(tokens[3]) + resx = float(tokens[4]) + resy = float(tokens[5]) + 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/gr3dnc01b.tif', cols, rows, 3, gdal.GDT_Float32) + ds.SetMetadataItem( + 'TIFFTAG_COPYRIGHT', 'Derived from work by Service Topographique, DITTT, Government of New Caledonia. Open License https://www.etalab.gouv.fr/wp-content/uploads/2014/05/Open_Licence.pdf') + ds.SetMetadataItem('TIFFTAG_IMAGEDESCRIPTION', + 'Geocentric translation from IGN72 GRANDE TERRE (EPSG:4662) to RGNC91-93 (EPSG:4906). Converted from gr3dnc01b.mnt') + 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', 'New Caledonia') + ds.SetMetadataItem('source_crs_wkt', """GEODCRS["IGN72 Grande Terre", + DATUM["IGN72 Grande Terre", + ELLIPSOID["International 1924",6378388,297, + LENGTHUNIT["metre",1]], + ID["EPSG",6634]], + PRIMEM["Greenwich",0, + ANGLEUNIT["degree",0.0174532925199433], + ID["EPSG",8901]], + 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]]]""") + ds.SetMetadataItem('target_crs_epsg_code', '4906') + 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(4749) + ds.SetProjection(sr.ExportToWkt()) + else: + tokens = [] + for t in l.split(' '): + if t: + tokens.append(t) + assert len(tokens) == 6 + lon = float(tokens[0]) + lat = float(tokens[1]) + dx = float(tokens[2]) + dy = float(tokens[3]) + dz = float(tokens[4]) + iy = (line_no - 1) % rows + ix = (line_no - 1) // 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 + + line_no += 1 + + assert count_val_rows == rows * cols + gdal.GetDriverByName('GTiff').CreateCopy('gr3dnc01b.tif', ds, options=[ + 'COMPRESS=DEFLATE', 'PREDICTOR=3', 'INTERLEAVE=BAND']) diff --git a/nc_dittt/convert_gr3dnc02b.py b/nc_dittt/convert_gr3dnc02b.py new file mode 100755 index 0000000..bb4ed13 --- /dev/null +++ b/nc_dittt/convert_gr3dnc02b.py @@ -0,0 +1,138 @@ +#!/usr/bin/env python +############################################################################### +# $Id$ +# +# Project: PROJ +# Purpose: Convert New Caledonia geocentric gr3dnc02b.mnt grid +# Author: Even Rouault +# +############################################################################### +# Copyright (c) 2019, Even Rouault +# +# 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 + +with open('gr3dnc02b.mnt', 'rt') as f: + count_val_rows = 0 + line_no = 0 + for l in f.readlines(): + if l[-1] == '\n': + l = l[0:-1] + if l[-1] == '\r': + l = l[0:-1] + if not l: + continue + if line_no == 0: + tokens = [] + for t in l.split(' '): + if t: + tokens.append(t) + assert len(tokens) >= 13 + minx = float(tokens[0]) + maxx = float(tokens[1]) + miny = float(tokens[2]) + maxy = float(tokens[3]) + resx = float(tokens[4]) + resy = float(tokens[5]) + 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/gr3dnc02b.tif', cols, rows, 3, gdal.GDT_Float32) + ds.SetMetadataItem( + 'TIFFTAG_COPYRIGHT', 'Derived from work by Service Topographique, DITTT, Government of New Caledonia. Open License https://www.etalab.gouv.fr/wp-content/uploads/2014/05/Open_Licence.pdf') + ds.SetMetadataItem('TIFFTAG_IMAGEDESCRIPTION', + 'Geocentric translation from IGN72 GRANDE TERRE (EPSG:4662) to RGNC91-93 (EPSG:4906). Converted from gr3dnc01b.mnt') + 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', 'New Caledonia') + ds.SetMetadataItem('source_crs_wkt', """GEODCRS["IGN72 Grande Terre", + DATUM["IGN72 Grande Terre", + ELLIPSOID["International 1924",6378388,297, + LENGTHUNIT["metre",1]], + ID["EPSG",6634]], + PRIMEM["Greenwich",0, + ANGLEUNIT["degree",0.0174532925199433], + ID["EPSG",8901]], + 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]]]""") + ds.SetMetadataItem('target_crs_epsg_code', '4906') + 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(4749) + ds.SetProjection(sr.ExportToWkt()) + else: + tokens = [] + for t in l.split('\t'): + if t: + tokens.append(t) + assert len(tokens) == 6 + lon = float(tokens[0]) + lat = float(tokens[1]) + dx = float(tokens[2]) + dy = float(tokens[3]) + dz = float(tokens[4]) + iy = (line_no - 1) % rows + ix = (line_no - 1) // 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 + + line_no += 1 + + assert count_val_rows == rows * cols + gdal.GetDriverByName('GTiff').CreateCopy('gr3dnc02b.tif', ds, options=[ + 'COMPRESS=DEFLATE', 'PREDICTOR=3', 'INTERLEAVE=BAND']) diff --git a/nc_dittt/convert_gr3dnc03a.py b/nc_dittt/convert_gr3dnc03a.py new file mode 100755 index 0000000..c8db270 --- /dev/null +++ b/nc_dittt/convert_gr3dnc03a.py @@ -0,0 +1,138 @@ +#!/usr/bin/env python +############################################################################### +# $Id$ +# +# Project: PROJ +# Purpose: Convert New Caledonia geocentric gr3dnc03a.mnt grid +# Author: Even Rouault +# +############################################################################### +# Copyright (c) 2019, Even Rouault +# +# 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 + +with open('gr3dnc03a.mnt', 'rt') as f: + count_val_rows = 0 + line_no = 0 + for l in f.readlines(): + if l[-1] == '\n': + l = l[0:-1] + if l[-1] == '\r': + l = l[0:-1] + if not l: + continue + if line_no == 0: + tokens = [] + for t in l.split(' '): + if t: + tokens.append(t) + assert len(tokens) >= 13 + minx = float(tokens[0]) + maxx = float(tokens[1]) + miny = float(tokens[2]) + maxy = float(tokens[3]) + resx = float(tokens[4]) + resy = float(tokens[5]) + 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/gr3dnc03a.tif', cols, rows, 3, gdal.GDT_Float32) + ds.SetMetadataItem( + 'TIFFTAG_COPYRIGHT', 'Derived from work by Service Topographique, DITTT, Government of New Caledonia. Open License https://www.etalab.gouv.fr/wp-content/uploads/2014/05/Open_Licence.pdf') + ds.SetMetadataItem('TIFFTAG_IMAGEDESCRIPTION', + 'Geocentric translation from NEA74 Noumea (EPSG:4644) to RGNC91-93 (EPSG:4906). Converted from gr3dnc03a.mnt') + 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', 'New Caledonia') + ds.SetMetadataItem('source_crs_wkt', """GEODCRS["NEA74 Noumea", + DATUM["NEA74 Noumea", + ELLIPSOID["International 1924",6378388,297, + LENGTHUNIT["metre",1]], + ID["EPSG",6644]], + PRIMEM["Greenwich",0, + ANGLEUNIT["degree",0.0174532925199433], + ID["EPSG",8901]], + 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]]]""") + ds.SetMetadataItem('target_crs_epsg_code', '4906') + 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(4749) + ds.SetProjection(sr.ExportToWkt()) + else: + tokens = [] + for t in l.split(' '): + if t: + tokens.append(t) + assert len(tokens) == 6 + lon = float(tokens[0]) + lat = float(tokens[1]) + dx = float(tokens[2]) + dy = float(tokens[3]) + dz = float(tokens[4]) + iy = (line_no - 1) % rows + ix = (line_no - 1) // 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 + + line_no += 1 + + assert count_val_rows == rows * cols + gdal.GetDriverByName('GTiff').CreateCopy('gr3dnc03a.tif', ds, options=[ + 'COMPRESS=DEFLATE', 'PREDICTOR=3', 'INTERLEAVE=BAND']) diff --git a/nc_dittt/nc_dittt_README.txt b/nc_dittt/nc_dittt_README.txt new file mode 100644 index 0000000..e3a0809 --- /dev/null +++ b/nc_dittt/nc_dittt_README.txt @@ -0,0 +1,33 @@ +# nc_dittt_README.txt + +The files in this section result from the conversion of datasets originating +from [Gouvernement de la Nouvelle-Calédonie – DITTT](https://dittt.gouv.nc/circe-pour-la-nouvelle-caledonie) + +## Included grids + +### New Caledonia: IGN72 GRANDE TERRE (EPSG:4662) to RGNC91-93 (EPSG:4906) (using Geocentric correction) + +*Source*: Gouvernement de la Nouvelle-Calédonie – DITTT](https://dittt.gouv.nc/circe-pour-la-nouvelle-caledonie) +*Format*: GeoTIFF converted from .mnt file with convert_gr3dnc01b.py +*License*: [License Ouverte (in French)](https://www.etalab.gouv.fr/wp-content/uploads/2017/04/ETALAB-Licence-Ouverte-v2.0.pdf) / [Open License (English translation)](https://www.etalab.gouv.fr/wp-content/uploads/2014/05/Open_Licence.pdf) +*Used by*: [EPSG:9329 Transformation](https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:coordinateOperation:EPSG::9329) + +* nc_dittt_gr3dnc01b.tif + +### New Caledonia - Noumea: IGN72 GRANDE TERRE (EPSG:4662) to RGNC91-93 (EPSG:4906) (using Geocentric correction) + +*Source*: Gouvernement de la Nouvelle-Calédonie – DITTT](https://dittt.gouv.nc/circe-pour-la-nouvelle-caledonie) +*Format*: GeoTIFF converted from .mnt file with convert_gr3dnc02b.py +*License*: [License Ouverte (in French)](https://www.etalab.gouv.fr/wp-content/uploads/2017/04/ETALAB-Licence-Ouverte-v2.0.pdf) / [Open License (English translation)](https://www.etalab.gouv.fr/wp-content/uploads/2014/05/Open_Licence.pdf) +*Used by*: [EPSG:9330 Transformation](https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:coordinateOperation:EPSG::9330) + +* nc_dittt_gr3dnc02b.tif + +### New Caledonia - Noumea: NEA74 Noumea (EPSG:4644) to RGNC91-93 (EPSG:4906) (using Geocentric correction) + +*Source*: Gouvernement de la Nouvelle-Calédonie – DITTT](https://dittt.gouv.nc/circe-pour-la-nouvelle-caledonie) +*Format*: GeoTIFF converted from .mnt file with convert_gr3dnc03a.py +*License*: [License Ouverte (in French)](https://www.etalab.gouv.fr/wp-content/uploads/2017/04/ETALAB-Licence-Ouverte-v2.0.pdf) / [Open License (English translation)](https://www.etalab.gouv.fr/wp-content/uploads/2014/05/Open_Licence.pdf) +*Used by*: [EPSG:9328 Transformation](https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:coordinateOperation:EPSG::9328) + +* nc_dittt_gr3dnc03a.tif diff --git a/nc_dittt/nc_dittt_gr3dnc01b.tif b/nc_dittt/nc_dittt_gr3dnc01b.tif new file mode 100644 index 0000000..ca8c779 Binary files /dev/null and b/nc_dittt/nc_dittt_gr3dnc01b.tif differ diff --git a/nc_dittt/nc_dittt_gr3dnc02b.tif b/nc_dittt/nc_dittt_gr3dnc02b.tif new file mode 100644 index 0000000..96c2c89 Binary files /dev/null and b/nc_dittt/nc_dittt_gr3dnc02b.tif differ diff --git a/nc_dittt/nc_dittt_gr3dnc03a.tif b/nc_dittt/nc_dittt_gr3dnc03a.tif new file mode 100644 index 0000000..786da45 Binary files /dev/null and b/nc_dittt/nc_dittt_gr3dnc03a.tif differ diff --git a/travis/expected_main.lst b/travis/expected_main.lst index 801849b..5a39017 100644 --- a/travis/expected_main.lst +++ b/travis/expected_main.lst @@ -136,6 +136,10 @@ is_lmi_ISN2004_ISN2016.tif is_lmi_ISN93_ISN2016.tif is_lmi_ISN_vel_beta.tif is_lmi_README.txt +nc_dittt_gr3dnc01b.tif +nc_dittt_gr3dnc02b.tif +nc_dittt_gr3dnc03a.tif +nc_dittt_README.txt NKG nl_nsgi_naptrans2018.tif nl_nsgi_nlgeo2018.tif -- cgit v1.2.3