120 lines
4.0 KiB
Python
120 lines
4.0 KiB
Python
#!/usr/bin/env python3
|
|
###############################################################################
|
|
# $Id: tolatlong.py e4fe7cc06270e5f38dfe78e6785a6bcca4e39e29 2021-04-01 21:02:04 +0300 Idan Miara $
|
|
#
|
|
# Project: GDAL Python samples
|
|
# Purpose: Script to read coordinate system and geotransformation matrix
|
|
# from input file and report latitude/longitude coordinates for the
|
|
# specified pixel.
|
|
# Author: Andrey Kiselev, dron@remotesensing.org
|
|
#
|
|
###############################################################################
|
|
# Copyright (c) 2003, Andrey Kiselev <dron@remotesensing.org>
|
|
# Copyright (c) 2009-2010, Even Rouault <even dot 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.
|
|
###############################################################################
|
|
|
|
import sys
|
|
|
|
from osgeo import gdal
|
|
from osgeo import osr
|
|
|
|
# =============================================================================
|
|
|
|
|
|
def Usage():
|
|
print('')
|
|
print('Read coordinate system and geotransformation matrix from input')
|
|
print('file and report latitude/longitude coordinates for the center')
|
|
print('of the specified pixel.')
|
|
print('')
|
|
print('Usage: tolatlong.py pixel line infile')
|
|
print('')
|
|
return 1
|
|
|
|
|
|
def main(argv):
|
|
infile = None
|
|
pixel = None
|
|
line = None
|
|
|
|
# =============================================================================
|
|
# Parse command line arguments.
|
|
# =============================================================================
|
|
i = 1
|
|
while i < len(argv):
|
|
arg = argv[i]
|
|
|
|
if pixel is None:
|
|
pixel = float(arg)
|
|
|
|
elif line is None:
|
|
line = float(arg)
|
|
|
|
elif infile is None:
|
|
infile = arg
|
|
|
|
else:
|
|
return Usage()
|
|
|
|
i = i + 1
|
|
|
|
if infile is None:
|
|
return Usage()
|
|
if pixel is None:
|
|
return Usage()
|
|
if line is None:
|
|
return Usage()
|
|
|
|
# Open input dataset
|
|
indataset = gdal.Open(infile, gdal.GA_ReadOnly)
|
|
|
|
# Read geotransform matrix and calculate ground coordinates
|
|
geomatrix = indataset.GetGeoTransform()
|
|
X = geomatrix[0] + geomatrix[1] * pixel + geomatrix[2] * line
|
|
Y = geomatrix[3] + geomatrix[4] * pixel + geomatrix[5] * line
|
|
|
|
# Shift to the center of the pixel
|
|
X += geomatrix[1] / 2.0
|
|
Y += geomatrix[5] / 2.0
|
|
|
|
# Build Spatial Reference object based on coordinate system, fetched from the
|
|
# opened dataset
|
|
srs = osr.SpatialReference()
|
|
if srs.ImportFromWkt(indataset.GetProjection()) != 0:
|
|
print("ERROR: Cannot import projection '%s'" % indataset.GetProjection())
|
|
return 1
|
|
|
|
srsLatLong = srs.CloneGeogCS()
|
|
ct = osr.CoordinateTransformation(srs, srsLatLong)
|
|
(lon, lat, height) = ct.TransformPoint(X, Y)
|
|
|
|
# Report results
|
|
print('pixel: %g\t\t\tline: %g' % (pixel, line))
|
|
print('latitude: %fd\t\tlongitude: %fd' % (lat, lon))
|
|
print('latitude: %s\t\tlongitude: %s' % (gdal.DecToDMS(lat, 'Lat', 2), gdal.DecToDMS(lon, 'Long', 2)))
|
|
|
|
return 0
|
|
|
|
|
|
if __name__ == '__main__':
|
|
sys.exit(main(sys.argv))
|
|
|