Python osgeo.gdal.GDT_Byte() Examples

The following are 27 code examples of osgeo.gdal.GDT_Byte(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module osgeo.gdal , or try the search function .
Example #1
Source Project: lidar   Author: giswqs   File: slicing.py    License: MIT License 6 votes vote down vote up
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 
Example #2
Source Project: apls   Author: CosmiQ   File: apls_utils.py    License: Apache License 2.0 6 votes vote down vote up
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


############################################################################### 
Example #3
Source Project: TF-SegNet   Author: mathildor   File: areacoverimage.py    License: MIT License 6 votes vote down vote up
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 
Example #4
Source Project: DsgTools   Author: dsgoficial   File: raster_processing.py    License: GNU General Public License v2.0 6 votes vote down vote up
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 
Example #5
Source Project: DsgTools   Author: dsgoficial   File: HSV_fusion.py    License: GNU General Public License v2.0 6 votes vote down vote up
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 
Example #6
Source Project: FloodMapsWorkshop   Author: vightel   File: landsat8_toa_watermap.py    License: Apache License 2.0 6 votes vote down vote up
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 
Example #7
Source Project: FloodMapsWorkshop   Author: vightel   File: landsat8_composite_toa.py    License: Apache License 2.0 6 votes vote down vote up
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 
Example #8
Source Project: buzzard   Author: airware   File: _gdal_gdt_conv.py    License: Apache License 2.0 5 votes vote down vote up
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] 
Example #9
Source Project: lidar   Author: giswqs   File: slicing.py    License: MIT License 5 votes vote down vote up
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 
Example #10
Source Project: lidar   Author: giswqs   File: filling.py    License: MIT License 5 votes vote down vote up
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 
Example #11
Source Project: coded   Author: bullocke   File: classify.py    License: MIT License 5 votes vote down vote up
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 
Example #12
Source Project: computing-pipeline   Author: terraref   File: Get_FLIR.py    License: BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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)) 
Example #13
Source Project: DsgTools   Author: dsgoficial   File: raster_processing.py    License: GNU General Public License v2.0 5 votes vote down vote up
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 
Example #14
Source Project: DsgTools   Author: dsgoficial   File: raster_processing.py    License: GNU General Public License v2.0 5 votes vote down vote up
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()) 
Example #15
Source Project: DsgTools   Author: dsgoficial   File: HSV_fusion.py    License: GNU General Public License v2.0 5 votes vote down vote up
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 
Example #16
Source Project: FloodMapsWorkshop   Author: vightel   File: viirs.py    License: Apache License 2.0 5 votes vote down vote up
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) 
Example #17
Source Project: FloodMapsWorkshop   Author: vightel   File: digiglobe_rgb_watermap.py    License: Apache License 2.0 5 votes vote down vote up
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 
Example #18
Source Project: FloodMapsWorkshop   Author: vightel   File: digiglobe_rgb_watermap.py    License: Apache License 2.0 5 votes vote down vote up
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 
Example #19
Source Project: FloodMapsWorkshop   Author: vightel   File: landsat8_toa_watermap.py    License: Apache License 2.0 5 votes vote down vote up
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 
Example #20
Source Project: FloodMapsWorkshop   Author: vightel   File: landsat8_composite_toa.py    License: Apache License 2.0 5 votes vote down vote up
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 
Example #21
Source Project: FloodMapsWorkshop   Author: vightel   File: palsar2.py    License: Apache License 2.0 4 votes vote down vote up
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 
Example #22
Source Project: FloodMapsWorkshop   Author: vightel   File: frost_pgc.py    License: Apache License 2.0 4 votes vote down vote 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) 
Example #23
Source Project: FloodMapsWorkshop   Author: vightel   File: hand_floodfill.py    License: Apache License 2.0 4 votes vote down vote up
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 
Example #24
Source Project: FloodMapsWorkshop   Author: vightel   File: modis_browseimage.py    License: Apache License 2.0 4 votes vote down vote up
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 
Example #25
Source Project: FloodMapsWorkshop   Author: vightel   File: eo1_to_topojson.py    License: Apache License 2.0 4 votes vote down vote up
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
	# 
Example #26
Source Project: FloodMapsWorkshop   Author: vightel   File: landsat8_composite_toa.py    License: Apache License 2.0 4 votes vote down vote up
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 
Example #27
Source Project: BlenderGIS   Author: domlysz   File: bigtiffwriter.py    License: GNU General Public License v3.0 4 votes vote down vote up
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')