PostGIS  2.5.7dev-r@@SVN_REVISION@@
ovdump.py
Go to the documentation of this file.
1 #! /usr/bin/env python
2 #
3 #
4 # This is a simple script based on GDAL to dump overview to separate file.
5 # It is used in WKTRaster testing to compare raster samples.
6 #
7 # NOTE: GDAL does include Frank's (unofficial) dumpoverviews utility too,
8 # which dumps overview as complete geospatial raster dataset.
9 #
10 # Copyright (C) 2009 Mateusz Loskot <mateusz@loskot.net>
11 #
12 # This program is free software; you can redistribute it and/or
13 # modify it under the terms of the GNU General Public License
14 # as published by the Free Software Foundation; either version 2
15 # of the License, or (at your option) any later version.
16 #
17 # This program is distributed in the hope that it will be useful,
18 # but WITHOUT ANY WARRANTY; without even the implied warranty of
19 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 # GNU General Public License for more details.
21 #
22 # You should have received a copy of the GNU General Public License
23 # along with this program; if not, write to the Free Software Foundation,
24 # Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
25 #
26 
27 from osgeo import gdal
28 from osgeo import osr
29 import osgeo.gdalconst as gdalc
30 from optparse import OptionParser
31 import os
32 import struct
33 import sys
34 
35 
36 prs = OptionParser(version="%prog $Revision: 4037 $",
37  usage="%prog -r <RASTER> -v <OVERVIEW>",
38  description="Dump GDAL raster overview to separate file, GeoTIF by default")
39 prs.add_option("-r", "--raster", dest="raster", action="store", default=None,
40  help="input raster file")
41 prs.add_option("-v", "--overview", dest="overview", action="store", type="int", default=None,
42  help="1-based index of raster overview, optional")
43 
44 (opts, args) = prs.parse_args()
45 
46 if opts.raster is None:
47  prs.error("use option -r to specify input raster file that contains overviews")
48 
49 
51 
52 # Source
53 src_ds = gdal.Open(opts.raster, gdalc.GA_ReadOnly);
54 if src_ds is None:
55  sys.exit('ERROR: Cannot open input file: ' + opts.raster)
56 
57 band = src_ds.GetRasterBand(1)
58 if opts.overview is None:
59  nvstart = 0
60  nvend = band.GetOverviewCount()
61 else:
62  nvstart = opts.overview - 1
63  nvend = opts.overview
64 
65 for nv in range(nvstart, nvend):
66 
67  band = src_ds.GetRasterBand(1)
68  if nv < 0 and nv >= band.GetOverviewCount():
69  print "ERROR: Failed to find overview %d" % nv
70  sys.exit(1)
71  ov_band = band.GetOverview(nv)
72 
73  ovf = int(0.5 + band.XSize / float(ov_band.XSize))
74 
75  print "--- OVERVIEW #%d level = %d ---------------------------------------" % (nv + 1, ovf)
76 
77  # Destination
78  dst_file = os.path.basename(opts.raster) + "_ov_" + str(ovf) + ".tif"
79  dst_format = "GTiff"
80  dst_drv = gdal.GetDriverByName(dst_format)
81  dst_ds = dst_drv.Create(dst_file, ov_band.XSize, ov_band.YSize, src_ds.RasterCount, ov_band.DataType)
82 
83  print "Source: " + opts.raster
84  print "Target: " + dst_file
85  print "Exporting %d bands of size %d x %d" % (src_ds.RasterCount, ov_band.XSize, ov_band.YSize)
86 
87  # Rewrite bands of overview to output file
88  ov_band = None
89  band = None
90 
91  for nb in range(1, src_ds.RasterCount + 1):
92 
93  band = src_ds.GetRasterBand(nb)
94  assert band is not None
95  ov_band = band.GetOverview(nv)
96  assert ov_band is not None
97 
98  print " Band #%d (%s) (%d x %d)" % \
99  (nb, gdal.GetDataTypeName(ov_band.DataType), ov_band.XSize, ov_band.YSize)
100 
101  dst_band = dst_ds.GetRasterBand(nb)
102  assert dst_band is not None
103  data = ov_band.ReadRaster(0, 0, ov_band.XSize, ov_band.YSize)
104  assert data is not None
105  dst_band.WriteRaster(0, 0, ov_band.XSize, ov_band.YSize, data, buf_type = ov_band.DataType)
106 
107  dst_ds = None
108 
109 # Cleanup
110 src_ds = None