Python osgeo.gdal.GDT_Float32() Examples

The following are 25 code examples of osgeo.gdal.GDT_Float32(). 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: buzzard   Author: airware   File: tools.py    License: Apache License 2.0 6 votes vote down vote up
def make_tif(path, tloffset=(0, 0), reso=(0.25, -0.25), rsize=(20, 10),
             proj=SRS[0]['wkt'], channel_count=1, dtype=gdal.GDT_Float32):
    """Create a tiff files and return info about it"""
    tl = ROOT_TL + tloffset
    reso = np.asarray(reso)
    fp = buzz.Footprint(tl=tl, rsize=rsize, size=np.abs(reso * rsize))
    x, y = fp.meshgrid_spatial
    x = np.abs(x) - abs(ROOT_TL[0])
    y = abs(ROOT_TL[1]) - np.abs(y)
    x *= 15
    y *= 15
    a = x / 2 + y / 2
    a = np.around(a).astype('float32')
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(path, rsize[0], rsize[1], channel_count, dtype)
    dataset.SetGeoTransform(fp.gt)
    dataset.SetProjection(proj)
    for i in range(channel_count):
        dataset.GetRasterBand(i + 1).WriteArray(a)
        dataset.GetRasterBand(i + 1).SetNoDataValue(-32000.)
    dataset.FlushCache()
    return path, fp, a 
Example #2
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 #3
Source Project: wradlib   Author: wradlib   File: test_io.py    License: MIT License 6 votes vote down vote up
def test_gdal_create_dataset(self):
        testfunc = io.gdal.gdal_create_dataset
        tmp = tempfile.NamedTemporaryFile(mode="w+b").name
        with pytest.raises(TypeError):
            testfunc("AIG", tmp)
        from osgeo import gdal

        with pytest.raises(TypeError):
            testfunc(
                "AAIGrid", tmp, cols=10, rows=10, bands=1, gdal_type=gdal.GDT_Float32
            )
        testfunc("GTiff", tmp, cols=10, rows=10, bands=1, gdal_type=gdal.GDT_Float32)
        testfunc(
            "GTiff",
            tmp,
            cols=10,
            rows=10,
            bands=1,
            gdal_type=gdal.GDT_Float32,
            remove=True,
        ) 
Example #4
Source Project: pyMeteo   Author: Mo-Dabao   File: tools.py    License: GNU Affero General Public License v3.0 6 votes vote down vote up
def write_geotiff(tiff_path, data, upper_left_lon, upper_left_lat, step, AREA_OR_POINT='Point'):
    """
    写入geotiff。默认pixel is point,默认WGS84
    """
    rows, columns = data.shape
    bands = 1
    pixel_width = step
    pixel_height = -step
    half_step = step / 2
    upper_left_x = upper_left_lon - half_step
    upper_left_y = upper_left_lat + half_step
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(tiff_path, columns, rows, bands, gdal.GDT_Float32)
    dataset.SetMetadataItem('AREA_OR_POINT', AREA_OR_POINT)
    dataset.SetGeoTransform([upper_left_x, pixel_width, 0,
                             upper_left_y, 0, pixel_height])
    # 获取地理坐标系统信息,用于选取需要的地理坐标系统
    sr = osr.SpatialReference()
    sr.SetWellKnownGeogCS('WGS84')
    dataset.SetProjection(sr.ExportToWkt())  # 给新建图层赋予投影信息
    dataset.GetRasterBand(1).WriteArray(data)
    dataset.FlushCache()  # 将数据写入硬盘
    dataset = None 
Example #5
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 #6
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 #7
Source Project: georaster   Author: GeoUtils   File: georaster.py    License: MIT License 6 votes vote down vote up
def save_geotiff(self,filename, dtype=gdal.GDT_Float32, **kwargs):
        """
        Save a GeoTIFF of the raster currently loaded.

        Georeferenced and subset according to the current raster.

        :params filename: the path and filename to save the file to.
        :type filename: str
        :params dtype: GDAL datatype, defaults to Float32.
        :type dtype: int

        :returns: 1 on success.
        :rtype: int

        """

        if self.r is None:
            raise ValueError('There are no image data loaded. No file can be created.')

        simple_write_geotiff(filename, self.r, self.trans, 
            wkt=self.srs.ExportToWkt(), dtype=dtype, **kwargs) 
Example #8
Source Project: pysteps   Author: pySTEPS   File: exporters.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def _create_geotiff_file(outfn, driver, shape, metadata, num_bands):
    dst = driver.Create(
        outfn,
        shape[1],
        shape[0],
        num_bands,
        gdal.GDT_Float32,
        ["COMPRESS=DEFLATE", "PREDICTOR=3"],
    )

    sx = (metadata["x2"] - metadata["x1"]) / shape[1]
    sy = (metadata["y2"] - metadata["y1"]) / shape[0]
    dst.SetGeoTransform([metadata["x1"], sx, 0.0, metadata["y2"], 0.0, -sy])

    sr = osr.SpatialReference()
    sr.ImportFromProj4(metadata["projection"])
    dst.SetProjection(sr.ExportToWkt())

    return dst 
Example #9
Source Project: buzzard   Author: airware   File: tools.py    License: Apache License 2.0 5 votes vote down vote up
def make_tif2(path, reso=(1., -1.), rsize=(10, 10), tl=(0., 10.),
              proj=SRS[0]['wkt'], channel_count=1, dtype=gdal.GDT_Float32,
              nodata=-32000, nodata_border_size=(0, 0, 0, 0)):
    """Create a tiff files"""
    reso = np.asarray(reso)
    fp = buzz.Footprint(tl=tl, rsize=rsize, size=np.abs(reso * rsize))
    x, y = fp.meshgrid_raster
    a = x + y
    if nodata_border_size != 0:
        l, r, t, b = nodata_border_size
        if t != 0:
            a[None:t, None:None] = nodata
        if b != 0:
            a[-b:None, None:None] = nodata
        if l != 0:
            a[None:None, None:l] = nodata
        if r != 0:
            a[None:None, -r:None] = nodata

    LOGGER.info('TIFF ARRAY:%s\n', a)
    gdal.UseExceptions()
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(path, int(rsize[0]), int(rsize[1]), channel_count, dtype)
    dataset.SetGeoTransform(fp.gt)
    dataset.SetProjection(proj)
    for i in range(channel_count):
        dataset.GetRasterBand(i + 1).WriteArray(a)
        dataset.GetRasterBand(i + 1).SetNoDataValue(nodata)
    dataset.FlushCache() 
Example #10
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 #11
Source Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_GDALIO.py    License: MIT License 5 votes vote down vote up
def array2raster(rasterfn,newRasterfn,array,driver_name = "ENVI", noDataValue = -9999):
    """Takes an array and writes to a GDAL compatible raster. It needs another raster to map the dimensions.

    Args:
        FileName (str): The filename (with path and extension) of a raster that has the same dimensions as the raster to be written.
        newRasterfn (str): The filename (with path and extension) of the new raster.
        array (np.array): The array to be written
        driver_name (str): The type of raster to write. Default is ENVI since that is the LSDTOpoTools format
        noDataValue (float): The no data value

    Return:
        np.array: A numpy array with the data from the raster.

    Author: SMM
    """

    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    cols = raster.RasterXSize
    rows = raster.RasterYSize

    driver = gdal.GetDriverByName(driver_name)
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outRaster.GetRasterBand(1).SetNoDataValue( noDataValue )
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()
#============================================================================== 
Example #12
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 #13
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 #14
Source Project: unmixing   Author: arthur-e   File: utils.py    License: MIT License 5 votes vote down vote up
def dump_raster(rast, rast_path, driver='GTiff', gdt=None, nodata=None):
    '''
    Creates a raster file from a given gdal.Dataset instance. Arguments:
        rast        A gdal.Dataset; does NOT accept NumPy array
        rast_path   The path of the output raster file
        driver      The name of the GDAL driver to use (determines file type)
        gdt         The GDAL data type to use, e.g., see gdal.GDT_Float32
        nodata      The NoData value; defaults to -9999.
    '''
    if gdt is None:
        gdt = rast.GetRasterBand(1).DataType
    driver = gdal.GetDriverByName(driver)
    sink = driver.Create(
        rast_path, rast.RasterXSize, rast.RasterYSize, rast.RasterCount, int(gdt))
    assert sink is not None, 'Cannot create dataset; there may be a problem with the output path you specified'
    sink.SetGeoTransform(rast.GetGeoTransform())
    sink.SetProjection(rast.GetProjection())

    for b in range(1, rast.RasterCount + 1):
        dat = rast.GetRasterBand(b).ReadAsArray()
        sink.GetRasterBand(b).WriteArray(dat)
        sink.GetRasterBand(b).SetStatistics(*map(np.float64,
            [dat.min(), dat.max(), dat.mean(), dat.std()]))

        if nodata is None:
            nodata = rast.GetRasterBand(b).GetNoDataValue()

            if nodata is None:
                nodata = -9999

        sink.GetRasterBand(b).SetNoDataValue(np.float64(nodata))

    sink.FlushCache() 
Example #15
Source Project: scikit-image-clustering-scripts   Author: h4k1m0u   File: io.py    License: MIT License 5 votes vote down vote up
def write(arr, filename):
        """
        Write image array to a file with the given filename.

        Args:
            arr (numpy.ndarray)
            filename (str)
        """
        driver = gdal.GetDriverByName('GTiff')
        x_size = arr.shape[1]
        y_size = arr.shape[0]
        dataset = driver.Create(filename, x_size, y_size,
                                eType=gdal.GDT_Float32)
        dataset.GetRasterBand(1).WriteArray(arr) 
Example #16
Source Project: FloodMapsWorkshop   Author: vightel   File: eo1_ali_l1t.py    License: Apache License 2.0 5 votes vote down vote up
def write_data(self, data, fileName, type=gdal.GDT_Float32, nbands=1, ct=0):
		fileName 	= os.path.join(self.outpath, fileName)
		
		if self.verbose:
			print "write_data", fileName
		
		driver 		= gdal.GetDriverByName( "GTiff" )
		dst_ds 		= driver.Create( fileName, self.RasterXSize, self.RasterYSize, nbands, type, [ 'INTERLEAVE=PIXEL', 'COMPRESS=DEFLATE' ] )
		band 		= dst_ds.GetRasterBand(1)

		band.WriteArray(data, 0, 0)
		band.SetNoDataValue(0)

		if ct :
			if self.verbose:
				print "write ct"
			
			ct = self.generate_color_table()
			band.SetRasterColorTable(ct)
			
		if self.verbose:
			print "write geotransform and projection"
		dst_ds.SetGeoTransform( self.geotransform )
		dst_ds.SetProjection( self.projection )
			
		if self.verbose:
			print "Written", fileName

		dst_ds 		= None 
Example #17
Source Project: RHEAS   Author: nasa   File: dbio.py    License: MIT License 5 votes vote down vote up
def writeGeotif(lat, lon, res, data, filename=None):
    """Writes Geotif in temporary directory so it can be imported into the PostGIS database."""
    if isinstance(data, np.ma.masked_array):
        nodata = np.double(data.fill_value)
        data = data.data
    else:
        nodata = -9999.
    if len(data.shape) < 2:
        nrows = int((max(lat) - min(lat)) / res) + 1
        ncols = int((max(lon) - min(lon)) / res) + 1
        out = np.zeros((nrows, ncols)) + nodata
        for c in range(len(lat)):
            i = int((max(lat) - lat[c]) / res)
            j = int((lon[c] - min(lon)) / res)
            out[i, j] = data[c]
    else:
        nrows, ncols = data.shape
        out = data
    if filename is None:
        f = tempfile.NamedTemporaryFile(suffix=".tif", delete=False)
        filename = f.name
        f.close()
    driver = gdal.GetDriverByName("GTiff")
    ods = driver.Create(filename, ncols, nrows, 1, gdal.GDT_Float32)
    ods.SetGeoTransform([min(lon) - res / 2.0, res, 0,
                         max(lat) + res / 2.0, 0, -res])
    srs = osr.SpatialReference()
    srs.SetWellKnownGeogCS("WGS84")
    ods.SetProjection(srs.ExportToWkt())
    ods.GetRasterBand(1).WriteArray(out)
    ods.GetRasterBand(1).SetNoDataValue(nodata)
    ods = None
    return filename 
Example #18
Source Project: RHEAS   Author: nasa   File: vic.py    License: MIT License 5 votes vote down vote up
def _writeRaster(self, data, filename):
        """Writes GeoTIFF raster temporarily so that it can be imported into the database."""
        nrows, ncols = data.shape
        driver = gdal.GetDriverByName("GTiff")
        ods = driver.Create(filename, ncols, nrows, 1, gdal.GDT_Float32)
        ods.SetGeoTransform([min(self.lon) - self.res / 2.0, self.res,
                             0, max(self.lat) + self.res / 2.0, 0, -self.res])
        srs = osr.SpatialReference()
        srs.SetWellKnownGeogCS("WGS84")
        ods.SetProjection(srs.ExportToWkt())
        ods.GetRasterBand(1).WriteArray(data)
        ods.GetRasterBand(1).SetNoDataValue(self.nodata)
        ods = None 
Example #19
Source Project: PyRate   Author: GeoscienceAustralia   File: shared.py    License: Apache License 2.0 5 votes vote down vote up
def gdal_dataset(out_fname, columns, rows, driver="GTiff", bands=1,
                 dtype='float32', metadata=None, crs=None,
                 geotransform=None, creation_opts=None):
    """
    Initialises a py-GDAL dataset object for writing image data.
    """
    if dtype == 'float32':
        gdal_dtype = gdal.GDT_Float32
    elif dtype == 'int16':
        gdal_dtype = gdal.GDT_Int16
    else:
        # assume gdal.GDT val is passed to function
        gdal_dtype = dtype

    # create output dataset
    driver = gdal.GetDriverByName(driver)
    outds = driver.Create(out_fname, columns, rows, bands, gdal_dtype, options=creation_opts)

    # geospatial info
    outds.SetGeoTransform(geotransform)
    outds.SetProjection(crs)

    # add metadata
    if metadata is not None:
        for k, v in metadata.items():
            outds.SetMetadataItem(k, str(v))

    return outds 
Example #20
Source Project: PyRate   Author: GeoscienceAustralia   File: shared.py    License: Apache License 2.0 5 votes vote down vote up
def write_output_geotiff(md, gt, wkt, data, dest, nodata):
    # pylint: disable=too-many-arguments
    """
    Writes PyRate output data to a GeoTIFF file.

    :param dict md: Dictionary containing PyRate metadata
    :param list gt: GDAL geotransform for the data
    :param list wkt: GDAL projection information for the data
    :param ndarray data: Output data array to save
    :param str dest: Destination file name
    :param float nodata: No data value of data

    :return None, file saved to disk
    """

    driver = gdal.GetDriverByName("GTiff")
    nrows, ncols = data.shape
    ds = driver.Create(dest, ncols, nrows, 1, gdal.GDT_Float32, options=['compress=packbits'])
    # set spatial reference for geotiff
    ds.SetGeoTransform(gt)
    ds.SetProjection(wkt)
    ds.SetMetadataItem(ifc.EPOCH_DATE, str(md[ifc.EPOCH_DATE]))

    # set other metadata
    ds.SetMetadataItem('DATA_TYPE', str(md['DATA_TYPE']))

    # sequence position for time series products

    for k in [ifc.SEQUENCE_POSITION, ifc.PYRATE_REFPIX_X, ifc.PYRATE_REFPIX_Y, ifc.PYRATE_REFPIX_LAT,
              ifc.PYRATE_REFPIX_LON, ifc.PYRATE_MEAN_REF_AREA, ifc.PYRATE_STDDEV_REF_AREA]:
        if k in md:
            ds.SetMetadataItem(k, str(md[k]))

    # write data to geotiff
    band = ds.GetRasterBand(1)
    band.SetNoDataValue(nodata)
    band.WriteArray(data, 0, 0) 
Example #21
Source Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_GDALIO.py    License: MIT License 4 votes vote down vote up
def RasterDifference(RasterFile1, RasterFile2, raster_band=1, OutFileName="Test.outfile", OutFileType="ENVI"):
    """
    Takes two rasters of same size and subtracts second from first,
    e.g. Raster1 - Raster2 = raster_of_difference
    then writes it out to file
    """

    Raster1 = gdal.Open(RasterFile1)
    Raster2 = gdal.Open(RasterFile2)

    print("RASTER 1: ")
    print(Raster1.GetGeoTransform())
    print(Raster1.RasterCount)
    print(Raster1.GetRasterBand(1).XSize)
    print(Raster1.GetRasterBand(1).YSize)
    print(Raster1.GetRasterBand(1).DataType)

    print("RASTER 2: ")
    print(Raster2.GetGeoTransform())
    print(Raster2.RasterCount)
    print(Raster2.GetRasterBand(1).XSize)
    print(Raster2.GetRasterBand(1).YSize)
    print(Raster2.GetRasterBand(1).DataType)

    raster_array1 = np.array(Raster1.GetRasterBand(raster_band).ReadAsArray())
    raster_array2 = np.array(Raster2.GetRasterBand(raster_band).ReadAsArray())

    assert(raster_array1.shape == raster_array2.shape )
    print("Shapes: ", raster_array1.shape, raster_array2.shape)

    difference_raster_array = raster_array1 - raster_array2

#    import matplotlib.pyplot as plt
#
#    plt.imshow(difference_raster_array)
#

    driver = gdal.GetDriverByName(OutFileType)

    dsOut = driver.Create(OutFileName,
                          Raster1.GetRasterBand(1).XSize,
                          Raster1.GetRasterBand(1).YSize,
                          1,
                          gdal.GDT_Float32)
                          #Raster1.GetRasterBand(raster_band).DataType)


    gdal_array.CopyDatasetInfo(Raster1,dsOut)
    bandOut = dsOut.GetRasterBand(1)
    gdal_array.BandWriteArray(bandOut, difference_raster_array)

#============================================================================== 
Example #22
Source Project: wradlib   Author: wradlib   File: gdal.py    License: MIT License 4 votes vote down vote up
def gdal_create_dataset(
    drv, name, cols=0, rows=0, bands=0, gdal_type=gdal.GDT_Unknown, remove=False
):
    """Creates GDAL.DataSet object.

    Parameters
    ----------
    drv : string
        GDAL driver string
    name : string
        path to filename
    cols : int
        # of columns
    rows : int
        # of rows
    bands : int
        # of raster bands
    gdal_type : raster data type
        eg. gdal.GDT_Float32
    remove : bool
        if True, existing gdal.Dataset will be
        removed before creation

    Returns
    -------
    out : gdal.Dataset
        object

    """
    driver = gdal.GetDriverByName(drv)
    metadata = driver.GetMetadata()

    if not metadata.get("DCAP_CREATE", False):
        raise TypeError(
            "WRADLIB: Driver {} doesn't support " "Create() method.".format(drv)
        )

    if remove:
        if os.path.exists(name):
            driver.Delete(name)
    ds = driver.Create(name, cols, rows, bands, gdal_type)

    return ds 
Example #23
Source Project: wa   Author: wateraccounting   File: raster_conversions.py    License: Apache License 2.0 4 votes vote down vote up
def Vector_to_Raster(Dir, shapefile_name, reference_raster_data_name):
    """
    This function creates a raster of a shp file

    Keyword arguments:
    Dir --
        str: path to the basin folder
    shapefile_name -- 'C:/....../.shp'
        str: Path from the shape file
    reference_raster_data_name -- 'C:/....../.tif'
        str: Path to an example tiff file (all arrays will be reprojected to this example)
    """

    from osgeo import gdal, ogr

    geo, proj, size_X, size_Y=Open_array_info(reference_raster_data_name)

    x_min = geo[0]
    x_max = geo[0] + size_X * geo[1]
    y_min = geo[3] + size_Y * geo[5]
    y_max = geo[3]
    pixel_size = geo[1]

    # Filename of the raster Tiff that will be created
    Dir_Basin_Shape = os.path.join(Dir,'Basin')
    if not os.path.exists(Dir_Basin_Shape):
        os.mkdir(Dir_Basin_Shape)

    Basename = os.path.basename(shapefile_name)
    Dir_Raster_end = os.path.join(Dir_Basin_Shape, os.path.splitext(Basename)[0]+'.tif')

    # Open the data source and read in the extent
    source_ds = ogr.Open(shapefile_name)
    source_layer = source_ds.GetLayer()

    # Create the destination data source
    x_res = int(round((x_max - x_min) / pixel_size))
    y_res = int(round((y_max - y_min) / pixel_size))

    # Create tiff file
    target_ds = gdal.GetDriverByName('GTiff').Create(Dir_Raster_end, x_res, y_res, 1, gdal.GDT_Float32, ['COMPRESS=LZW'])
    target_ds.SetGeoTransform(geo)
    srse = osr.SpatialReference()
    srse.SetWellKnownGeogCS(proj)
    target_ds.SetProjection(srse.ExportToWkt())
    band = target_ds.GetRasterBand(1)
    target_ds.GetRasterBand(1).SetNoDataValue(-9999)
    band.Fill(-9999)

    # Rasterize the shape and save it as band in tiff file
    gdal.RasterizeLayer(target_ds, [1], source_layer, None, None, [1], ['ALL_TOUCHED=TRUE'])
    target_ds = None

    # Open array
    Raster_Basin = Open_tiff_array(Dir_Raster_end)

    return(Raster_Basin) 
Example #24
Source Project: georaster   Author: GeoUtils   File: georaster.py    License: MIT License 4 votes vote down vote up
def from_array(cls, raster, geo_transform, proj4, 
        gdal_dtype=gdal.GDT_Float32, nodata=None):
        """ Create a georaster object from numpy array and georeferencing information.

        :param raster: 2-D NumPy array of raster to load
        :type raster: np.array
        :param geo_transform: a Geographic Transform tuple of the form \
        (top left x, w-e cell size, 0, top left y, 0, n-s cell size (-ve))
        :type geo_transform: tuple
        :param proj4: a proj4 string representing the raster projection
        :type proj4: str
        :param gdal_dtype: a GDAL data type (default gdal.GDT_Float32)
        :type gdal_dtype: int
        :param nodata: None or the nodata value for this array
        :type nodata: None, int, float, str

        :returns: GeoRaster object
        :rtype: GeoRaster

         """

        if len(raster.shape) > 2:
            nbands = raster.shape[2]
        else:
            nbands = 1

        # Create a GDAL memory raster to hold the input array
        mem_drv = gdal.GetDriverByName('MEM')
        source_ds = mem_drv.Create('', raster.shape[1], raster.shape[0],
                         nbands, gdal_dtype)

        # Set geo-referencing
        source_ds.SetGeoTransform(geo_transform)
        srs = osr.SpatialReference()
        srs.ImportFromProj4(proj4)
        source_ds.SetProjection(srs.ExportToWkt())

        # Write input array to the GDAL memory raster
        for b in range(0,nbands):
            if nbands > 1:
                r = raster[:,:,b]
            else:
                r = raster
            source_ds.GetRasterBand(b+1).WriteArray(r)
            if nodata != None:
                source_ds.GetRasterBand(b+1).SetNoDataValue(nodata)

        # Return a georaster instance instantiated by the GDAL raster
        return cls(source_ds) 
Example #25
Source Project: pysteps   Author: pySTEPS   File: exporters.py    License: BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
def _export_geotiff(F, exporter):
    def init_band(band):
        band.SetScale(1.0)
        band.SetOffset(0.0)
        band.SetUnitType(exporter["metadata"]["unit"])

    if exporter["incremental"] is None:
        for i in range(exporter["num_timesteps"]):
            if exporter["num_ens_members"] == 1:
                band = exporter["dst"][i].GetRasterBand(1)
                init_band(band)
                band.WriteArray(F[i, :, :])
            else:
                for j in range(exporter["num_ens_members"]):
                    band = exporter["dst"][i].GetRasterBand(j + 1)
                    init_band(band)
                    band.WriteArray(F[j, i, :, :])
    elif exporter["incremental"] == "timestep":
        i = exporter["num_files_written"]

        outfn = _get_geotiff_filename(
            exporter["outfnprefix"],
            exporter["startdate"],
            exporter["num_timesteps"],
            exporter["timestep"],
            i,
        )
        dst = _create_geotiff_file(
            outfn,
            exporter["driver"],
            exporter["shape"],
            exporter["metadata"],
            exporter["num_ens_members"],
        )

        for j in range(exporter["num_ens_members"]):
            band = dst.GetRasterBand(j + 1)
            init_band(band)
            if exporter["num_ens_members"] > 1:
                band.WriteArray(F[j, :, :])
            else:
                band.WriteArray(F)

        exporter["num_files_written"] += 1
    elif exporter["incremental"] == "member":
        for i in range(exporter["num_timesteps"]):
            # NOTE: This does not work because the GeoTIFF driver does not
            # support adding bands. An alternative solution needs to be
            # implemented.
            exporter["dst"][i].AddBand(gdal.GDT_Float32)
            band = exporter["dst"][i].GetRasterBand(exporter["dst"][i].RasterCount)
            init_band(band)
            band.WriteArray(F[i, :, :])