Python osgeo.gdal.GDT_Byte() Examples

def writeRaster(arr, out_path, template):
    no_data = 0
    # First of all, gather some information from the template file
    data = gdal.Open(template)
    [cols, rows] = arr.shape
    trans = data.GetGeoTransform()
    proj = data.GetProjection()
    # nodatav = 0 #data.GetNoDataValue()
    # Create the file, using the information from the template file
    outdriver = gdal.GetDriverByName("GTiff")
    # http://www.gdal.org/gdal_8h.html
    # GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6,
    outdata   = outdriver.Create(str(out_path), rows, cols, 1, gdal.GDT_UInt32)
    # Write the array to the file, which is the original array in this example
    outdata.GetRasterBand(1).WriteArray(arr)
    # Set a no data value if required
    outdata.GetRasterBand(1).SetNoDataValue(no_data)
    # Georeference the image
    outdata.SetGeoTransform(trans)
    # Write projection information
    outdata.SetProjection(proj)
    return arr


# raster to vector 

def write_data(self, data, fileName, colorTable):
		fileName 	= os.path.join(self.outpath, fileName)
		driver 		= gdal.GetDriverByName( "GTiff" )
		dst_ds 		= driver.Create( fileName, self.RasterXSize, self.RasterYSize, 1, gdal.GDT_Byte, [ 'COMPRESS=DEFLATE' ] )
		band 		= dst_ds.GetRasterBand(1)

		if self.geotransform:
			dst_ds.SetGeoTransform( self.geotransform )
			
		if self.projection:
			dst_ds.SetProjection( self.projection )

		if colorTable:
			print "Add colortable"
			band.SetRasterColorTable(self.ct)

		band.WriteArray(data, 0, 0)
			
		if verbose:
			print "Written", fileName

		ds 		= None 

def write_data(self, data, fileName, colorTable):
		fileName 	= os.path.join(self.outpath, fileName)
		driver 		= gdal.GetDriverByName( "GTiff" )
		dst_ds 		= driver.Create( fileName, self.RasterXSize, self.RasterYSize, 1, gdal.GDT_Byte, [ 'COMPRESS=DEFLATE' ] )
		band 		= dst_ds.GetRasterBand(1)

		if self.geotransform:
			dst_ds.SetGeoTransform( self.geotransform )
			
		if self.projection:
			dst_ds.SetProjection( self.projection )

		if colorTable:
			print "Add colortable"
			band.SetRasterColorTable(self.ct)

		band.WriteArray(data, 0, 0)
			
		if verbose:
			print "Written", fileName

		ds 		= None 

def getNumpyType(self, pixelType = gdal.GDT_Byte):
        """
        Translates the gdal raster type to numpy type
        pixelType: gdal raster type
        """
        if pixelType == gdal.GDT_Byte:
            return numpy.uint8
        elif pixelType == gdal.GDT_UInt16:
            return numpy.uint16
        elif pixelType == gdal.GDT_Int16:
            return numpy.int16
        elif pixelType == gdal.GDT_UInt32:
            return numpy.uint32
        elif pixelType == gdal.GDT_Int32:
            return numpy.int32
        elif pixelType == gdal.GDT_Float32:
            return numpy.float32
        elif pixelType == gdal.GDT_Float64:
            return numpy.float64 

def getNumpyType(self, pixelType = gdal.GDT_Byte):
        '''
        Translates the gdal raster type to numpy type
        pixelType: gdal raster type
        '''
        if pixelType == gdal.GDT_Byte:
            return numpy.uint8
        elif pixelType == gdal.GDT_UInt16:
            return numpy.uint16
        elif pixelType == gdal.GDT_Int16:
            return numpy.int16
        elif pixelType == gdal.GDT_UInt32:
            return numpy.uint32
        elif pixelType == gdal.GDT_Int32:
            return numpy.int32
        elif pixelType == gdal.GDT_Float32:
            return numpy.float32
        elif pixelType == gdal.GDT_Float64:
            return numpy.float64 

def createImageDS(filename, x_min, y_min, x_max, y_max, pixel_size,  srs=None):
    # Create the destination data source
    x_res = int((x_max - x_min) / pixel_size) # resolution
    y_res = int((y_max - y_min) / pixel_size) # resolution
    ds = gdal.GetDriverByName('GTiff').Create(filename, x_res,
            y_res, 1, gdal.GDT_Byte)
    ds.SetGeoTransform((
            x_min, pixel_size, 0,
            y_max, 0, -pixel_size,
        ))
    if srs:
        # Make the target raster have the same projection as the source
        ds.SetProjection(srs.ExportToWkt())
    else:
        # Source has no projection (needs GDAL >= 1.7.0 to work)
        ds.SetProjection('LOCAL_CS["arbitrary"]')

    # Set nodata
    band = ds.GetRasterBand(1)
    band.SetNoDataValue(0)

    return ds 

def CreateMultiBandGeoTiff(OutPath, Array):
    '''
    Author: Jake Shermeyer
    Array has shape:
        Channels, Y, X?
    '''
    driver = gdal.GetDriverByName('GTiff')
    DataSet = driver.Create(OutPath, Array.shape[2], Array.shape[1],
                            Array.shape[0], gdal.GDT_Byte,
                            ['COMPRESS=LZW'])
    for i, image in enumerate(Array, 1):
        DataSet.GetRasterBand(i).WriteArray(image)
    del DataSet

    return OutPath


############################################################################### 

def polygonize(img,shp_path):
    # mapping between gdal type and ogr field type
    type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                    gdal.GDT_UInt16: ogr.OFTInteger,
                    gdal.GDT_Int16: ogr.OFTInteger,
                    gdal.GDT_UInt32: ogr.OFTInteger,
                    gdal.GDT_Int32: ogr.OFTInteger,
                    gdal.GDT_Float32: ogr.OFTReal,
                    gdal.GDT_Float64: ogr.OFTReal,
                    gdal.GDT_CInt16: ogr.OFTInteger,
                    gdal.GDT_CInt32: ogr.OFTInteger,
                    gdal.GDT_CFloat32: ogr.OFTReal,
                    gdal.GDT_CFloat64: ogr.OFTReal}

    ds = gdal.Open(img)
    prj = ds.GetProjection()
    srcband = ds.GetRasterBand(1)
    dst_layername = "Shape"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference(wkt=prj)

    dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
    raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# convert images in a selected folder to shapefiles 

def create_geotiff(np_arr, gps_bounds, out_file_path, base_dir):
    try:
        nrows,ncols = np.shape(np_arr)
        # gps_bounds: (lat_min, lat_max, lng_min, lng_max)
        xres = (gps_bounds[3] - gps_bounds[2])/float(ncols)
        yres = (gps_bounds[1] - gps_bounds[0])/float(nrows)
        geotransform = (gps_bounds[2],xres,0,gps_bounds[1],0,-yres) #(top left x, w-e pixel resolution, rotation (0 if North is up), top left y, rotation (0 if North is up), n-s pixel resolution)

        output_path = out_file_path

        output_raster = gdal.GetDriverByName('GTiff').Create(output_path, ncols, nrows, 3, gdal.GDT_Byte)

        output_raster.SetGeoTransform(geotransform) # specify coordinates
        srs = osr.SpatialReference() # establish coordinate encoding
        srs.ImportFromEPSG(4326) # specifically, google mercator
        output_raster.SetProjection( srs.ExportToWkt() ) # export coordinate system to file

        # TODO: Something wonky w/ uint8s --> ending up w/ lots of gaps in data (white pixels)
        output_raster.GetRasterBand(1).WriteArray(np_arr.astype('uint8')) # write red channel to raster file
        output_raster.GetRasterBand(1).FlushCache()
        output_raster.GetRasterBand(1).SetNoDataValue(-99)
        
        output_raster.GetRasterBand(2).WriteArray(np_arr.astype('uint8')) # write green channel to raster file
        output_raster.GetRasterBand(2).FlushCache()
        output_raster.GetRasterBand(2).SetNoDataValue(-99)

        output_raster.GetRasterBand(3).WriteArray(np_arr.astype('uint8')) # write blue channel to raster file
        output_raster.GetRasterBand(3).FlushCache()
        output_raster.GetRasterBand(3).SetNoDataValue(-99)


        # for test: once we've saved the image, make sure to append this path to our list of TIFs
        tif_file = os.path.join(base_dir, tif_list_file)
        f = open(tif_file,'a+')
        f.write(output_path + '\n')
    except Exception as ex:
        fail('Error creating GeoTIFF: ' + str(ex)) 

def readBlock(self, band, sizeX, sizeY = 1, offsetY = 0, pixelType = gdal.GDT_Byte):
        '''
        Reads image block
        band: band used
        sizeX: x block size
        sizeY: y block size
        offsetY: Y offset
        pixelType: pixel type
        '''
        numpytype = self.getNumpyType(pixelType)

        bandscanline =band.ReadRaster( 0, offsetY, sizeX, sizeY, sizeX, sizeY, pixelType )
        pixelArray=numpy.fromstring(bandscanline, dtype=numpytype)
        
        return pixelArray 

def write_geotiff(fname, data, geo_transform, projection, data_type=gdal.GDT_Byte):

    """
    Create a GeoTIFF file with the given data.
    :param fname: Path to a directory with shapefiles
    :param data: Number of rows of the result
    :param geo_transform: Returned value of 
	gdal.Dataset.GetGeoTransform (coefficients for transforming between 
	pixel/line (P,L) raster space, and projection coordinates (Xp,Yp) space.
    :param projection: Projection definition string (Returned by 
	gdal.Dataset.GetProjectionRef)
    """

    driver = gdal.GetDriverByName('GTiff')
    rows, cols = data.shape
    dataset = driver.Create(fname, cols, rows, 1, data_type)
    dataset.SetGeoTransform(geo_transform)
    dataset.SetProjection(projection)
    band = dataset.GetRasterBand(1)
    band.WriteArray(data)

    metadata = {
        'TIFFTAG_COPYRIGHT': 'CC BY 4.0',
        'TIFFTAG_DOCUMENTNAME': 'classification',
        'TIFFTAG_IMAGEDESCRIPTION': 'Supervised classification.',
        'TIFFTAG_SOFTWARE': 'Python, GDAL, scikit-learn'
    }
    dataset.SetMetadata(metadata)

    dataset = None  # Close the file
    return 

def writeBlock(self, band, block, sizeX, sizeY = 1, offsetY = 0, pixelType = gdal.GDT_Byte):
        '''
        Writes image block
        band: band used
        sizeX: x block size
        sizeY: y block size
        offsetY: Y offset
        pixelType: pixel type
        '''
        numpytype = self.getNumpyType(pixelType)

        band.WriteRaster(0, offsetY, sizeX, sizeY, block.astype(numpytype).tostring()) 

def readBlock(self, band, sizeX, sizeY = 1, offsetY = 0, pixelType = gdal.GDT_Byte):
        """
        Reads image block
        band: band used
        sizeX: x block size
        sizeY: y block size
        offsetY: Y offset
        pixelType: pixel type
        """
        numpytype = self.getNumpyType(pixelType)

        bandscanline =band.ReadRaster( 0, offsetY, sizeX, sizeY, sizeX, sizeY, pixelType )
        pixelArray=numpy.fromstring(bandscanline, dtype=numpytype)
        
        return pixelArray 

def polygonize(img,shp_path):
    # mapping between gdal type and ogr field type
    type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                    gdal.GDT_UInt16: ogr.OFTInteger,
                    gdal.GDT_Int16: ogr.OFTInteger,
                    gdal.GDT_UInt32: ogr.OFTInteger,
                    gdal.GDT_Int32: ogr.OFTInteger,
                    gdal.GDT_Float32: ogr.OFTReal,
                    gdal.GDT_Float64: ogr.OFTReal,
                    gdal.GDT_CInt16: ogr.OFTInteger,
                    gdal.GDT_CInt32: ogr.OFTInteger,
                    gdal.GDT_CFloat32: ogr.OFTReal,
                    gdal.GDT_CFloat64: ogr.OFTReal}

    ds = gdal.Open(img)
    prj = ds.GetProjection()
    srcband = ds.GetRasterBand(1)

    dst_layername = "Shape"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference(wkt=prj)

    dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
    # raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType])
    raster_field = ogr.FieldDefn('id', type_mapping[gdal.GDT_Int32])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# extract sinks from dem 

def _str_of_gdt(gdt):
    return {
        gdal.GDT_Byte: 'GDT_Byte',
        gdal.GDT_Int16: 'GDT_Int16',
        gdal.GDT_Int32: 'GDT_Int32',
        gdal.GDT_UInt16: 'GDT_UInt16',
        gdal.GDT_UInt32: 'GDT_UInt32',
        gdal.GDT_Float32: 'GDT_Float32',
        gdal.GDT_Float64: 'GDT_Float64',
        gdal.GDT_CFloat32: 'GDT_CFloat32',
        gdal.GDT_CFloat64: 'GDT_CFloat64',
        gdal.GDT_CInt16: 'GDT_CInt16',
        gdal.GDT_CInt32: 'GDT_CInt32',
    }[gdt] 

def CreateTopojsonFile(self, fileName, src_ds, projection, geotransform, ct, data, pres, xorg, ymax, water ):
	
		driver 				= gdal.GetDriverByName( "GTiff" )
		dst_ds_dataset		= driver.Create( fileName, src_ds.RasterXSize, src_ds.RasterYSize, 1, gdal.GDT_Byte, [ 'COMPRESS=DEFLATE' ] )
		
		dst_ds_dataset.SetGeoTransform( geotransform )
		dst_ds_dataset.SetProjection( projection )

		o_band				= dst_ds_dataset.GetRasterBand(1)
	
		o_band.SetRasterColorTable(ct)
		o_band.WriteArray(data, 0, 0)

		dst_ds_dataset = None
		print "Created", fileName

		cmd = "gdal_translate -q -of PNM -expand gray " + fileName + " "+fileName+".bmp"
		self.execute(cmd)

		# -i  		invert before processing
		# -t 2  	suppress speckles of up to this many pixels. 
		# -a 1.5  	set the corner threshold parameter
		# -z black  specify how to resolve ambiguities in path decomposition. Must be one of black, white, right, left, minority, majority, or random. Default is minority
		# -x 		scaling factor
		# -L		left margin
		# -B		bottom margin

		cmd = str.format("potrace -i -z black -a 1.5 -t 3 -b geojson -o {0} {1} -x {2} -L {3} -B {4} ", fileName+".geojson", fileName+".bmp", pres, xorg, ymax ); 
		self.execute(cmd)
	
		cmd = str.format("topojson -o {0} --simplify-proportion 0.75 -p water={1} -- water={2}", fileName+".topojson", water, fileName+".geojson" ); 
		self.execute(cmd)
	
		# convert it back to json
		cmd = "topojson-geojson --precision 4 -o %s %s" % ( self.geojsonDir, fileName+".topojson" )
		self.execute(cmd)
	
		# rename file
		output_file = "water_level_%d.geojson" % water
		cmd = "mv %s %s" % (os.path.join(self.geojsonDir,"water.json"), os.path.join(self.geojsonDir, output_file))
		self.execute(cmd) 

def computeTOAReflectance(self, band, data):
		mp 			= float(self.metadata['REFLECTANCE_MULT_BAND_'+str(band)])
		ap			= float(self.metadata['REFLECTANCE_ADD_BAND_'+str(band)])
		se			= float(self.metadata['SUN_ELEVATION'])
		
		if verbose:
			print 'REFLECTANCE_MULT_BAND_'+str(band), mp
			print 'REFLECTANCE_ADD_BAND_'+str(band), ap
			print 'SUN_ELEVATION', se, math.sin( se * math.pi/180.0)
			
		toa			= (mp * data + ap) / math.sin( se * math.pi/180.0)
		
		toa[ toa < 0 ] = 0
		
		# save it for debug purpose
		if 0 :
			fname 		= os.path.join(self.outpath, "band_" + str(band) + "_toa.tif")
			driver 		= gdal.GetDriverByName( "GTiff" )
			dst_ds 		= driver.Create( fname, self.RasterXSize, self.RasterYSize, 1, gdal.GDT_Byte )
			band 		= dst_ds.GetRasterBand(1)			
			band.WriteArray(self.linear_stretch(toa), 0, 0)
			dst_ds.SetGeoTransform( self.geotransform )
			dst_ds.SetProjection( self.projection )
		
			dst_ds		= None
			print "Written TOA", fname, self.RasterXSize, self.RasterYSize, numpy.min(toa), numpy.max(toa), numpy.mean(toa), numpy.std(toa) 
		
		return toa 

def gray1( self, red, green, blue, RasterXSize, RasterYSize ):		
		lum			= 0.95 * red + 0.01 * green #+ 0.0 * blue
		driver 		= gdal.GetDriverByName( "GTiff" )
		dst_ds 		= driver.Create( self.gray_file, RasterXSize, RasterYSize, 1, gdal.GDT_Byte, [ 'INTERLEAVE=PIXEL', 'COMPRESS=DEFLATE' ] )
		band 		= dst_ds.GetRasterBand(1)
		band.WriteArray(lum, 0, 0)
		band.SetNoDataValue(0)
		dst_ds 		= None
		print "Written", self.gray_file
		return lum
	
	# http://en.wikipedia.org/wiki/Relative_luminance 

def rgb2gray( self, red, green, blue, RasterXSize, RasterYSize ):		
		lum			= 0.2126 * red + 0.7152 * green + 0.0722 * blue
		
		#driver 		= gdal.GetDriverByName( "GTiff" )
		#dst_ds 		= driver.Create( self.output_file2, RasterXSize, RasterYSize, 1, gdal.GDT_Byte, [ 'INTERLEAVE=PIXEL', 'COMPRESS=DEFLATE' ] )
		#band 		= dst_ds.GetRasterBand(1)
		#band.WriteArray(lum, 0, 0)
		#band.SetNoDataValue(0)
		#dst_ds 		= None
		#print "Written", self.output_file2
		return lum 

def computeTOAReflectance(self, band, data):
		mp 			= float(self.metadata['REFLECTANCE_MULT_BAND_'+str(band)])
		ap			= float(self.metadata['REFLECTANCE_ADD_BAND_'+str(band)])
		se			= float(self.metadata['SUN_ELEVATION'])
		
		if verbose:
			print 'REFLECTANCE_MULT_BAND_'+str(band), mp
			print 'REFLECTANCE_ADD_BAND_'+str(band), ap
			print 'SUN_ELEVATION', se, math.sin( se * math.pi/180.0)
			
		toa			= (mp * data + ap) / math.sin( se * math.pi/180.0)
		
		toa[ toa < 0 ] = 0
		
		# save it for debug purpose
		if 0 :
			fname 		= os.path.join(self.outpath, "band_" + str(band) + "_toa.tif")
			driver 		= gdal.GetDriverByName( "GTiff" )
			dst_ds 		= driver.Create( fname, self.RasterXSize, self.RasterYSize, 1, gdal.GDT_Byte )
			band 		= dst_ds.GetRasterBand(1)			
			band.WriteArray(self.linear_stretch(toa), 0, 0)
			dst_ds.SetGeoTransform( self.geotransform )
			dst_ds.SetProjection( self.projection )
		
			dst_ds		= None
			print "Written TOA", fname, self.RasterXSize, self.RasterYSize, numpy.min(toa), numpy.max(toa), numpy.mean(toa), numpy.std(toa) 
		
		return toa 

def process(self):
			
		red_data  	= self.get_band_data(self.red)
		red_data[self.cloud_mask]	=	0
		red_data[self.cirrus_mask]	=	0		
		red_data	= self.linear_stretch(red_data)
		
		green_data  = self.get_band_data(self.green)
		green_data[self.cloud_mask]	=	0
		green_data[self.cirrus_mask]=	0		
		green_data	= self.linear_stretch(green_data)
			
		blue_data  	= self.get_band_data(self.blue)
		blue_data[self.cloud_mask]	=	0
		blue_data[self.cirrus_mask]	=	0		
		blue_data	= self.linear_stretch(blue_data)
		
		if verbose:
			print "Creating", self.output_file
			
		driver 			= gdal.GetDriverByName( "GTiff" )
		dst_ds 			= driver.Create( self.output_file, self.RasterXSize, self.RasterYSize, 4, gdal.GDT_Byte, [ 'INTERLEAVE=PIXEL', 'COMPRESS=DEFLATE' ] )
		
		if verbose:
			print "Writing red band"
			
		dst_ds.GetRasterBand(1).WriteArray(red_data, 0, 0)
		dst_ds.GetRasterBand(1).SetNoDataValue(0)
		
		if verbose:
			print "Writing green band"
		dst_ds.GetRasterBand(2).WriteArray(green_data, 0, 0)
		dst_ds.GetRasterBand(2).SetNoDataValue(0)
		
		if verbose:
			print "Writing blue band"
		dst_ds.GetRasterBand(3).WriteArray(blue_data, 0, 0)
		dst_ds.GetRasterBand(3).SetNoDataValue(0)
		
		if verbose:
			print "Writing alpha band"
		alpha_band			= dst_ds.GetRasterBand(4)
		alpha_data 			= alpha_band.ReadAsArray(0, 0, dst_ds.RasterXSize, dst_ds.RasterYSize )
		alpha_data[self.no_data]	= 0
		alpha_data[blue_data>0]		= 255
		alpha_band.WriteArray(alpha_data, 0, 0)

		if self.geotransform:
			dst_ds.SetGeoTransform( self.geotransform )
			
		if self.projection:
			dst_ds.SetProjection( self.projection )

		dst_ds 	= None
		ds		= None 

def __init__(self, path, w, h, georef, geoTiffOptions={'TFW':'YES', 'TILED':'YES', 'BIGTIFF':'YES', 'COMPRESS':'JPEG', 'JPEG_QUALITY':80, 'PHOTOMETRIC':'YCBCR'}):
		'''
		path = fule system path for the ouput tiff
		w, h = width and height in pixels
		georef : a Georef object used to set georeferencing informations, optional
		geoTiffOptions : GDAL create option for tiff format
		'''

		if not HAS_GDAL:
			raise ImportError("GDAL interface unavailable")


		#control path validity

		self.w = w
		self.h = h
		self.size = (w, h)

		self.path = path
		self.georef = georef

		if geoTiffOptions.get('COMPRESS', None) == 'JPEG':
			#JPEG in tiff cannot have an alpha band, workaround is to use internal tiff mask
			self.useMask = True
			gdal.SetConfigOption('GDAL_TIFF_INTERNAL_MASK', 'YES')
			n = 3 #RGB
		else:
			self.useMask = False
			n = 4 #RGBA
		self.nbBands = n

		options = [str(k) + '=' + str(v) for k, v in geoTiffOptions.items()]

		driver = gdal.GetDriverByName("GTiff")
		gdtype = gdal.GDT_Byte #GDT_UInt16, GDT_Int16, GDT_UInt32, GDT_Int32
		self.dtype = 'uint8'

		self.ds = driver.Create(path, w, h, n, gdtype, options)
		if self.useMask:
			self.ds.CreateMaskBand(gdal.GMF_PER_DATASET)#The mask band is shared between all bands on the dataset
			self.mask = self.ds.GetRasterBand(1).GetMaskBand()
			self.mask.Fill(255)
		elif n == 4:
			self.ds.GetRasterBand(4).Fill(255)

		#Write georef infos
		self.ds.SetGeoTransform(self.georef.toGDAL())
		if self.georef.crs is not None:
			self.ds.SetProjection(self.georef.crs.getOgrSpatialRef().ExportToWkt())
		#self.georef.toWorldFile(os.path.splitext(path)[0] + '.tfw') 

def save_hand_as_png(self):
		hand_ds = gdal.Open( self.hand_file )
		if hand_ds is None:
			print('ERROR: hand file no data:')
			sys.exit(-1)

		RasterXSize = hand_ds.RasterXSize
		RasterYSize = hand_ds.RasterYSize
		RasterCount = hand_ds.RasterCount
		
		print "Hand file", RasterXSize, RasterYSize, RasterCount
		
		driver 		= gdal.GetDriverByName( "GTiff" )
		dst_ds 		= driver.Create( self.hand_rgb, RasterXSize, RasterYSize, 1, gdal.GDT_Byte,
			[ 'INTERLEAVE=PIXEL', 'COMPRESS=DEFLATE' ] )

		ct = gdal.ColorTable()
		for i in range(256):
			ct.SetColorEntry( i, (255, 255, 255, 255) )

		# Colorbrewer, sequential, 7
		ct.SetColorEntry( 0, (0, 0, 0, 255) )
		ct.SetColorEntry( 1, (8, 48, 107, 255) )
		ct.SetColorEntry( 2, (8, 81, 156, 255) )
		ct.SetColorEntry( 3, (33, 113, 181, 255) )
		ct.SetColorEntry( 4, (66, 146, 198, 255) )
		ct.SetColorEntry( 5, (107, 174, 214, 255) )
		ct.SetColorEntry( 6, (158, 202, 225, 255) )
		ct.SetColorEntry( 7, (198, 219, 239, 255) )
		ct.SetColorEntry( 8, (222, 235, 2247, 255) )
		ct.SetColorEntry( 9, (247, 251, 255, 255) )

		# ocean
		ct.SetColorEntry( 255, (0, 0, 0, 0) )

		hand = hand_ds.GetRasterBand(1)
		data = hand.ReadAsArray(0, 0, RasterXSize, RasterYSize )
		
		band = dst_ds.GetRasterBand(1)
		band.SetRasterColorTable(ct)
		band.WriteArray(data, 0, 0)
		band.SetNoDataValue(0)
		
		# copy projection
		projection   = hand_ds.GetProjection()
		geotransform = hand_ds.GetGeoTransform()

		dst_ds.SetGeoTransform( geotransform )
		dst_ds.SetProjection( projection )

		dst_ds 		= None
		hand_ds 	= None
		
	#
	# Generate Hand file
	# 

def process(infile, coastlines, outfile):
	ds = gdal.Open( infile )
		
	# Surface Water
	band = ds.GetRasterBand(1)
	data = band.ReadAsArray(0, 0, ds.RasterXSize, ds.RasterYSize )
	data = (data >= 2)*255
	
	coastlines_ds	= gdal.Open(coastlines)
	coastal_band 	= coastlines_ds.GetRasterBand(1)
	#coastal_data 	= coastal_band.ReadAsArray(0, 0, coastlines_ds.RasterXSize, coastlines_ds.RasterYSize )
	coastal_data 	= coastal_band.ReadAsArray(0, 0, ds.RasterXSize, ds.RasterYSize )
	
	# Coastal Masking
	mask = coastal_data>0
	data[mask]= 0
	
	# Step 1
	# extract surface water from MWP product
	#
	driver 		= gdal.GetDriverByName( "GTIFF" )
	dst_ds 		= driver.Create( outfile, ds.RasterXSize, ds.RasterYSize, 4, gdal.GDT_Byte, [ 'INTERLEAVE=PIXEL', 'COMPRESS=DEFLATE' ] )

	#ct = gdal.ColorTable()
	#for i in range(256):
	#	ct.SetColorEntry( i, (255, 255, 255, 255) )
			
	#ct.SetColorEntry( 0, (0, 0, 0, 0) )
	#ct.SetColorEntry( 1, (255, 0, 0, 0) )
	#ct.SetColorEntry( 2, (255, 0, 0, 0) )
	#ct.SetColorEntry( 3, (255, 0, 0, 255) )
		
	band = dst_ds.GetRasterBand(1)
	#band.SetRasterColorTable(ct)
	band.WriteArray(data, 0, 0)
	band.SetNoDataValue(0)

	alphaband = dst_ds.GetRasterBand(4)
	alphaband.WriteArray(data, 0, 0)

	geotransform = ds.GetGeoTransform()
	projection   = ds.GetProjection()
	
	dst_ds.SetGeoTransform( geotransform )
	dst_ds.SetProjection( projection )

	dst_ds 			= None
	ds 				= None
	coastlines_ds	= None
	
#	
# Generates a browseimage
# browseimage.py -y 2013 -d 205 -t 020E010S -p 2 -v 

def save_data(self):
		print str(datetime.now()), "Saving Hand data..."
		print "  water pix:", self.wp
		print "  processed:", self.count
		print "  total pixs:", self.RasterXSize * self.RasterYSize
		print "  Hand pixs:", numpy.count_nonzero(self.hand)
		
		hand_img 	= os.path.join(self.inpath, self.zone, self.tile + "_hand.tif")
		driver 		= gdal.GetDriverByName( "GTiff" )
		#dst_ds 		= driver.CreateCopy( hand_img, self.hds, 0, [ 'TILED=YES', 'COMPRESS=PACKBITS' ] )
		dst_ds 		= driver.Create( hand_img, self.RasterXSize, self.RasterYSize, 1, gdal.GDT_Byte,
			[ 'INTERLEAVE=PIXEL', 'COMPRESS=DEFLATE' ] )

		ct = gdal.ColorTable()
		for i in range(256):
			ct.SetColorEntry( i, (255, 255, 255, 255) )

		# Colorbrewer, sequential, 7
		ct.SetColorEntry( 0, (0, 0, 0, 255) )
		ct.SetColorEntry( 1, (8, 48, 107, 255) )
		ct.SetColorEntry( 2, (8, 81, 156, 255) )
		ct.SetColorEntry( 3, (33, 113, 181, 255) )
		ct.SetColorEntry( 4, (66, 146, 198, 255) )
		ct.SetColorEntry( 5, (107, 174, 214, 255) )
		ct.SetColorEntry( 6, (158, 202, 225, 255) )
		ct.SetColorEntry( 7, (198, 219, 239, 255) )
		ct.SetColorEntry( 8, (222, 235, 2247, 255) )
		ct.SetColorEntry( 9, (247, 251, 255, 255) )
		ct.SetColorEntry( 254, (255, 0, 0, 255) )

		# ocean
		ct.SetColorEntry( 255, (0, 0, 0, 0) )

		band = dst_ds.GetRasterBand(1)
		band.SetRasterColorTable(ct)
		band.WriteArray(self.hand, 0, 0)
		band.SetNoDataValue(0)
		
		# copy projection
		projection   = self.hds.GetProjection()
		geotransform = self.hds.GetGeoTransform()

		dst_ds.SetGeoTransform( geotransform )
		dst_ds.SetProjection( projection )

		dst_ds 		= None
		self.hds 	= None

	# One row up 

def CreateTopojsonFile(srcPath, fileName, src_ds, projection, geotransform, ct, data, pres, xorg, ymax, frost ):
	
	geojsonDir			= os.path.join(srcPath,"geojson")
	if not os.path.exists(geojsonDir):            
		os.makedirs(geojsonDir)

	driver 				= gdal.GetDriverByName( "GTiff" )
	dst_ds_dataset		= driver.Create( fileName, src_ds.RasterXSize, src_ds.RasterYSize, 1, gdal.GDT_Byte, [ 'COMPRESS=DEFLATE' ] )
	dst_ds_dataset.SetGeoTransform( geotransform )
	dst_ds_dataset.SetProjection( projection )

	o_band				= dst_ds_dataset.GetRasterBand(1)
	
	o_band.SetRasterColorTable(ct)
	o_band.WriteArray(data, 0, 0)

	dst_ds_dataset = None
	print "Created", fileName

	cmd = "gdal_translate -q -of PNM -expand gray " + fileName + " "+fileName+".pgm"
	execute(cmd)

	# -i  		invert before processing
	# -t 2  	suppress speckles of up to this many pixels. 
	# -a 1.5  	set the corner threshold parameter
	# -z black  specify how to resolve ambiguities in path decomposition. Must be one of black, white, right, left, minority, majority, or random. Default is minority
	# -x 		scaling factor
	# -L		left margin
	# -B		bottom margin

	cmd = str.format("potrace -z black -a 1.5 -t 3 -b geojson -o {0} {1} -x {2} -L {3} -B {4} ", fileName+".geojson", fileName+".pgm", pres, xorg, ymax ); 
	execute(cmd)

	#cmd = str.format("node set_geojson_property.js --file {0} --prop frost={1}", fileName+".geojson", frost)
	#execute(cmd)
	
	cmd = str.format("topojson -o {0} --simplify-proportion 0.75 -p frost={1} -- frost={2}", fileName+".topojson", frost, fileName+".geojson" ); 
	execute(cmd)
	
	# convert it back to json
	cmd = "topojson-geojson --precision 4 -o %s %s" % ( geojsonDir, fileName+".topojson" )
	execute(cmd)
	
	# rename file
	output_file = "frost_level_%d.geojson" % frost
	cmd = "mv %s %s" % (os.path.join(geojsonDir,"frost.json"), os.path.join(geojsonDir, output_file))
	execute(cmd) 

def process_raw_data(self):
		if self.force or not os.path.isfile(self.rgb_scene):
			format 				= "GTiff"
			driver 				= gdal.GetDriverByName( format )

			if verbose:
				print "Reading", self.hh_scene 
				
			hh_ds 				= gdal.Open( self.hh_scene )
			hh_band 			= hh_ds.GetRasterBand(1)
			
			self.RasterXSize 	= hh_ds.RasterXSize
			self.RasterYSize 	= hh_ds.RasterYSize

			self.data			= hh_band.ReadAsArray(0, 0, self.RasterXSize, self.RasterYSize )
			
			projection  		= hh_ds.GetProjection()
			geotransform		= hh_ds.GetGeoTransform()
			
			self.threshold()
			#self.speckle_filter('median', 3)
			
			self.generate_color_table()

			# will be open as writeable as well
			if verbose:
				print "Creating output...", self.surface_water
				
			dst_ds				= driver.Create( self.surface_water, self.RasterXSize, self.RasterYSize, 1, gdal.GDT_Byte, [ 'COMPRESS=DEFLATE' ]  )

			band = dst_ds.GetRasterBand(1)
			
			print "set color table"
			band.SetRasterColorTable(self.ct)

			print "write array"
			band.WriteArray(self.data, 0, 0)

			print "set nodata"
			band.SetNoDataValue(0)
		
			dst_ds.SetGeoTransform( geotransform )
			dst_ds.SetProjection( projection )

			dst_ds 	= None
			hh_ds	= None
			hv_ds	= None
			
			print "Done!"
		
# palsar2.py --scene ALOS2004073570-140620 -v