|
| 1 | +#!/usr/bin/env python |
| 2 | +#------------------------------------------------------------------------------ |
| 3 | +# |
| 4 | +# This tool clips data to mask. Based on the mask value, the code |
| 5 | +# performs folloging pixel operations: |
| 6 | +# mask no-data value (0x00) -> sets pixel to a given no-data value |
| 7 | +# mask data value (0xFF) -> copies pixel from the original image |
| 8 | +# |
| 9 | +# This tool extracts subset of the image specified by the row/column |
| 10 | +# offset of the upper-left corenr and row/column size of extracted |
| 11 | +# block. The tool takes care about preserving the geo-metadata. |
| 12 | +# |
| 13 | +# Project: Image Processing Tools |
| 14 | +# Authors: Martin Paces <[email protected]> |
| 15 | +# |
| 16 | +#------------------------------------------------------------------------------- |
| 17 | +# Copyright (C) 2013 EOX IT Services GmbH |
| 18 | +# |
| 19 | +# Permission is hereby granted, free of charge, to any person obtaining a copy |
| 20 | +# of this software and associated documentation files (the "Software"), to deal |
| 21 | +# in the Software without restriction, including without limitation the rights |
| 22 | +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell |
| 23 | +# copies of the Software, and to permit persons to whom the Software is |
| 24 | +# furnished to do so, subject to the following conditions: |
| 25 | +# |
| 26 | +# The above copyright notice and this permission notice shall be included in all |
| 27 | +# copies of this Software or works derived from this Software. |
| 28 | +# |
| 29 | +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |
| 30 | +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
| 31 | +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE |
| 32 | +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
| 33 | +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, |
| 34 | +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN |
| 35 | +# THE SOFTWARE. |
| 36 | +#------------------------------------------------------------------------------- |
| 37 | + |
| 38 | +import sys |
| 39 | +import os.path |
| 40 | +import img_block as ib |
| 41 | +import img_util as iu |
| 42 | +import numpy as np |
| 43 | + |
| 44 | +def clipToMask( bi , bm , mask_fg , nodata ) : |
| 45 | + |
| 46 | + if bi.non_equal_2d( bm ) : |
| 47 | + raise RuntimeError( "Equal blocks' extents required!" ) |
| 48 | + |
| 49 | + # prepare list of no-data values |
| 50 | + if len( nodata ) == 1 : |
| 51 | + nodata = [ nodata[0] for i in xrange(bi.data.shape[2]) ] |
| 52 | + |
| 53 | + fg = ( bm.data[:,:,0] == mask_fg ) |
| 54 | + bg = ( bm.data[:,:,0] != mask_fg ) |
| 55 | + |
| 56 | + for i in xrange( bi.data.shape[2] ) : |
| 57 | + bi.data[:,:,i] = bg * nodata[i] + fg * bi.data[:,:,i] |
| 58 | + |
| 59 | +#------------------------------------------------------------------------------ |
| 60 | + |
| 61 | +if __name__ == "__main__" : |
| 62 | + |
| 63 | + # TODO: to improve CLI |
| 64 | + |
| 65 | + EXENAME = os.path.basename( sys.argv[0] ) |
| 66 | + |
| 67 | + def error( message ) : |
| 68 | + print >>sys.stderr, "ERROR: %s: %s\n" %( EXENAME, message ) |
| 69 | + |
| 70 | + # block size |
| 71 | + bsx , bsy = 256, 256 |
| 72 | + |
| 73 | + # default format options |
| 74 | + |
| 75 | + FOPTS = ib.FormatOptions() |
| 76 | + FOPTS["TILED"] = "YES" |
| 77 | + FOPTS["BLOCKXSIZE"] = "256" |
| 78 | + FOPTS["BLOCKYSIZE"] = "256" |
| 79 | + FOPTS["COMPRESS"] = "DEFLATE" |
| 80 | + |
| 81 | + try: |
| 82 | + |
| 83 | + INPUT = sys.argv[1] |
| 84 | + MASK = sys.argv[2] |
| 85 | + OUTPUT = sys.argv[3] |
| 86 | + NODATA = sys.argv[4].split(",") |
| 87 | + MASKBG = 0x00 |
| 88 | + MASKFG = 0xFF |
| 89 | + |
| 90 | + #anything else treated as a format option |
| 91 | + for opt in sys.argv[5:] : |
| 92 | + FOPTS.setOption( opt ) |
| 93 | + |
| 94 | + except IndexError : |
| 95 | + |
| 96 | + error("Not enough input arguments!") |
| 97 | + sys.stderr.write("USAGE: %s <input image> <mask> <output TIF> <no data value or list>\n"%EXENAME) |
| 98 | + sys.stderr.write("EXAMPLE: %s input.tif mask.tif output.tif 255,255,255\n"%EXENAME) |
| 99 | + sys.stderr.write("EXAMPLE: %s input.tif mask.tif output.tif 0\n"%EXENAME) |
| 100 | + sys.exit(1) |
| 101 | + |
| 102 | + # open input image |
| 103 | + imi = ib.ImgFileIn( INPUT ) |
| 104 | + |
| 105 | + # convert no-data values to the image's data type |
| 106 | + NODATA = map( np.dtype(imi.dtype).type, NODATA ) |
| 107 | + |
| 108 | + if len(NODATA) < imi.sz : |
| 109 | + NODATA=[ NODATA[0] for i in xrange(imi.sz) ] |
| 110 | + |
| 111 | + # open input mask |
| 112 | + imm = ib.ImgFileIn( MASK ) |
| 113 | + |
| 114 | + # check mask properties |
| 115 | + |
| 116 | + if imm.sz > 1 : |
| 117 | + error("Multiband mask not supported!") |
| 118 | + sys.exit(1) |
| 119 | + |
| 120 | + if imm.dtype != 'uint8' : |
| 121 | + error("Unsupported mask data type '%s'!"%imi.dtype) |
| 122 | + sys.exit(1) |
| 123 | + |
| 124 | + if not imm.equal_2d( imi ) : |
| 125 | + error("Input mask and image must have the same pixel" |
| 126 | + " size! image: (%d x %d) mask: (%d x %d)" %( |
| 127 | + imi.sy , imi.sx, imm.sy , imm.sx ) ) |
| 128 | + sys.exit(1) |
| 129 | + |
| 130 | + # creation parameters |
| 131 | + prm = { |
| 132 | + 'path' : OUTPUT, |
| 133 | + 'nrow' : imi.sy, |
| 134 | + 'ncol' : imi.sx, |
| 135 | + 'nband' : imi.sz, |
| 136 | + 'dtype' : imi.dtype, |
| 137 | + 'options' : FOPTS.getOptions(), |
| 138 | + 'nodata' : NODATA, |
| 139 | + } |
| 140 | + |
| 141 | + #print prm |
| 142 | + |
| 143 | + # geocoding |
| 144 | + if imi.ds.GetProjection() : |
| 145 | + prm['proj'] = imi.ds.GetProjection() |
| 146 | + prm['geotrn'] = imi.ds.GetGeoTransform() |
| 147 | + elif imi.ds.GetGCPProjection() : |
| 148 | + prm['proj'] = imi.ds.GetGCPProjection() |
| 149 | + prm['gcps'] = imi.ds.GetGCPs() |
| 150 | + |
| 151 | + # open output image |
| 152 | + imo = ib.createGeoTIFF( **prm ) |
| 153 | + |
| 154 | + #-------------------------------------------------------------------------- |
| 155 | + |
| 156 | + # initialize progress printer |
| 157 | + prg = iu.Progress( (1+(imi.sy-1)/bsy)*(1+(imi.sx-1)/bsx) ) |
| 158 | + |
| 159 | + print "Clipping image by a mask ..." |
| 160 | + |
| 161 | + for ty in xrange( 1 + (imi.sy-1)/bsy ) : |
| 162 | + for tx in xrange( 1 + (imi.sx-1)/bsx ) : |
| 163 | + |
| 164 | + # extent of the tile |
| 165 | + ex_ti = imi & ib.ImgExtent( (bsx,bsy,imi.sz) , (tx*bsx,ty*bsy,0) ) |
| 166 | + ex_tm = imi & ib.ImgExtent( (bsx,bsy,1) , (tx*bsx,ty*bsy,0) ) |
| 167 | + |
| 168 | + # allocate input image block |
| 169 | + bi = ib.ImgBlock( imi.dtype , extent = ex_ti ) |
| 170 | + |
| 171 | + # allocate mask block |
| 172 | + bm = ib.ImgBlock( imm.dtype , extent = ex_tm ) |
| 173 | + |
| 174 | + # load image block |
| 175 | + imi.read( bi ) |
| 176 | + |
| 177 | + # load mask block |
| 178 | + imm.read( bm ) |
| 179 | + |
| 180 | + # clip image block to mask |
| 181 | + clipToMask( bi , bm , MASKFG , NODATA ) |
| 182 | + |
| 183 | + # save image block |
| 184 | + imo.write( bi ) |
| 185 | + |
| 186 | + # print progress |
| 187 | + sys.stdout.write(prg.istr(1)) ; sys.stdout.flush() |
| 188 | + |
| 189 | + sys.stdout.write("\n") ; sys.stdout.flush() |
0 commit comments