← Back to branch summary

~ubuntu-branches/debian/sid/gdal/sid

« back to all changes in this revision

Viewing changes to swig/python/samples/val_at_coord.py

  • Committer: Package Import Robot
  • Author(s): Francesco Paolo Lovergine
  • Date: 2012-05-07 15:04:42 UTC
  • mfrom: (5.5.16 experimental)
  • Revision ID: package-import@ubuntu.com-20120507150442-2eks97loeh6rq005
Tags: 1.9.0-1
* Ready for sid, starting transition.
* All symfiles updated to latest builds.
* Added dh_numpy call in debian/rules to depend on numpy ABI.
* Policy bumped to 3.9.3, no changes required.

Show diffs side-by-side

added added

removed removed

Lines of Context:
swig/python/samples/val_at_coord.py
 
1
#!/usr/bin/env python
 
2
###############################################################################
 
3
# $Id: val_at_coord.py 20674 2010-09-22 17:50:33Z rouault $
 
4
#
 
5
# Project:  GDAL Python samples
 
6
# Purpose:  Outputs the value of the raster bands at a given
 
7
#           (longitude, latitude) or (X, Y) location.
 
8
# Author:   Even Rouault
 
9
#
 
10
###############################################################################
 
11
# Copyright (c) 2010, Even Rouault <even dot rouault at mines dash paris dot org>
 
12
# 
 
13
# Permission is hereby granted, free of charge, to any person obtaining a
 
14
# copy of this software and associated documentation files (the "Software"),
 
15
# to deal in the Software without restriction, including without limitation
 
16
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
 
17
# and/or sell copies of the Software, and to permit persons to whom the
 
18
# Software is furnished to do so, subject to the following conditions:
 
19
#
 
20
# The above copyright notice and this permission notice shall be included
 
21
# in all copies or substantial portions of the Software.
 
22
#
 
23
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 
24
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 
25
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 
26
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 
27
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 
28
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 
29
# DEALINGS IN THE SOFTWARE.
 
30
###############################################################################
 
31
 
 
32
try:
 
33
    from osgeo import gdal
 
34
    from osgeo import osr
 
35
except ImportError:
 
36
    import gdal
 
37
 
 
38
import sys
 
39
 
 
40
# =============================================================================
 
41
def Usage():
 
42
    print('Usage: val_at_coord.py [-display_xy] [longitude latitude | -coordtype=georef X Y] filename')
 
43
    print('')
 
44
    print('By default, the 2 first arguments are supposed to be the location')
 
45
    print('in longitude, latitude order. If -coordtype=georef is specified before')
 
46
    print('the next 2 values will be interpretated as the X and Y coordinates')
 
47
    print('in the dataset spatial reference system.')
 
48
    sys.exit( 1 )
 
49
 
 
50
# =============================================================================
 
51
 
 
52
display_xy = False
 
53
coordtype_georef = False
 
54
longitude = None
 
55
latitude = None
 
56
filename = None
 
57
 
 
58
# =============================================================================
 
59
# Parse command line arguments.
 
60
# =============================================================================
 
61
i = 1
 
62
while i < len(sys.argv):
 
63
    arg = sys.argv[i]
 
64
 
 
65
    if arg == '-coordtype=georef':
 
66
        coordtype_georef = True
 
67
 
 
68
    elif arg == '-display_xy':
 
69
        display_xy = True
 
70
 
 
71
    elif longitude is None:
 
72
        longitude = float(arg)
 
73
 
 
74
    elif latitude is None:
 
75
        latitude = float(arg)
 
76
 
 
77
    elif filename is None:
 
78
        filename = arg
 
79
 
 
80
    else:
 
81
        Usage()
 
82
 
 
83
    i = i + 1
 
84
 
 
85
if longitude is None:
 
86
    Usage()
 
87
if latitude is None:
 
88
    Usage()
 
89
if filename is None:
 
90
    filename()
 
91
 
 
92
# Open input dataset
 
93
ds = gdal.Open( filename, gdal.GA_ReadOnly )
 
94
if ds is None:
 
95
    print('Cannot open %s' % filename)
 
96
    sys.exit(1)
 
97
 
 
98
# Build Spatial Reference object based on coordinate system, fetched from the
 
99
# opened dataset
 
100
if coordtype_georef:
 
101
    X = longitude
 
102
    Y = latitude
 
103
else:
 
104
    srs = osr.SpatialReference()
 
105
    srs.ImportFromWkt(ds.GetProjection())
 
106
 
 
107
    srsLatLong = srs.CloneGeogCS()
 
108
    # Convert from (longitude,latitude) to projected coordinates
 
109
    ct = osr.CoordinateTransformation(srsLatLong, srs)
 
110
    (X, Y, height) = ct.TransformPoint(longitude, latitude)
 
111
 
 
112
# Read geotransform matrix and calculate corresponding pixel coordinates
 
113
geomatrix = ds.GetGeoTransform()
 
114
(success, inv_geometrix) = gdal.InvGeoTransform(geomatrix)
 
115
x = int(inv_geometrix[0] + inv_geometrix[1] * X + inv_geometrix[2] * Y)
 
116
y = int(inv_geometrix[3] + inv_geometrix[4] * X + inv_geometrix[5] * Y)
 
117
 
 
118
if display_xy:
 
119
    print('x=%d, y=%d' % (x, y))
 
120
 
 
121
if x < 0 or x >= ds.RasterXSize or y < 0 or y >= ds.RasterYSize:
 
122
    print('Passed coordinates are not in dataset extent')
 
123
    sys.exit(1)
 
124
 
 
125
res = ds.ReadAsArray(x,y,1,1)
 
126
if len(res.shape) == 2:
 
127
    print(res[0][0])
 
128
else:
 
129
    for val in res:
 
130
        print(val[0][0])