Python osgeo.gdal.GDT_Float32() Examples

The following are code examples for showing how to use osgeo.gdal.GDT_Float32(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: GeoPy   Author: aerler   File: simple_regrid.py    GNU General Public License v3.0 6 votes vote down vote up
def getProj(self, bands, dtype='float32', size=None):
    '''
    generic function that returns a gdal dataset, ready for use 
    '''
    # determine GDAL data type
    if dtype == 'float32': gdt = gdal.GDT_Float32
    # determine size
    if not size: size = self.size # should be default  
    # create GDAL dataset 
    dset = ramdrv.Create('', int(size[0]), int(size[1]), int(bands), int(gdt)) 
    #if bands > 6: # add more bands, if necessary
      #for i in xrange(bands-6): dset.AddBand()
    # N.B.: for some reason a dataset is always initialized with 6 bands
    # set projection parameters
    dset.SetGeoTransform(self.geotransform) # does the order matter?
    dset.SetProjection(self.projection.ExportToWkt()) # is .ExportToWkt() necessary?
    # return dataset
    return dset

## simple lat/lon geo-referencing system 
Example 2
Project: buzzard   Author: airware   File: tools.py    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 3
Project: lidar   Author: giswqs   File: slicing.py    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 4
Project: lidar   Author: giswqs   File: slicing.py    MIT License 6 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 5
Project: watools   Author: wateraccounting   File: data_conversions.py    Apache License 2.0 6 votes vote down vote up
def Save_as_MEM(data='', geo='', projection=''):
    """
    This function save the array as a memory file

    Keyword arguments:
    data -- [array], dataset of the geotiff
    geo -- [minimum lon, pixelsize, rotation, maximum lat, rotation,
            pixelsize], (geospatial dataset)
    projection -- interger, the EPSG code
    """
    # save as a geotiff
    driver = gdal.GetDriverByName("MEM")
    dst_ds = driver.Create('', int(data.shape[1]), int(data.shape[0]), 1,
                           gdal.GDT_Float32)
    srse = osr.SpatialReference()
    if projection == '':
        srse.SetWellKnownGeogCS("WGS84")
    else:
        srse.SetWellKnownGeogCS(projection)
    dst_ds.SetProjection(srse.ExportToWkt())
    dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
    dst_ds.SetGeoTransform(geo)
    dst_ds.GetRasterBand(1).WriteArray(data)
    return(dst_ds) 
Example 6
Project: surfclass   Author: Kortforsyningen   File: rasterio.py    MIT License 6 votes vote down vote up
def gdaltype_to_creationoptions(gdaltype):
    """Gets GDAL creation options suitable for a GDAL datatype.

    Args:
        gdaltype (int): GDAL datatype value. For instance `osgeo.gdal.GDT_Byte`

    Raises:
        NotImplementedError: If given GDAL data type is not supported.

    Returns:
        list of str: List of GDAL creation options.

    """
    if gdaltype in [
        gdal.GDT_Byte,
        gdal.GDT_Int16,
        gdal.GDT_Int32,
        gdal.GDT_UInt16,
        gdal.GDT_UInt32,
    ]:
        return gdal_int_options
    if gdaltype in [gdal.GDT_Float32, gdal.GDT_Float64]:
        return gdal_float_options
    raise NotImplementedError() 
Example 7
Project: surfclass   Author: Kortforsyningen   File: test_cli_prepare.py    MIT License 6 votes vote down vote up
def test_cli_prepare_extractfeatures(cli_runner, amplituderaster_filepath, tmp_path):
    args = f"prepare extractfeatures -b 727000 6171000 728000 6172000 -f mean -f var -n 5 -c reflect {amplituderaster_filepath} {tmp_path}"

    result = cli_runner.invoke(cli, args.split(" "), catch_exceptions=False)
    assert result.exit_code == 0

    outfile = tmp_path / "mean.tif"
    ds = gdal.Open(str(outfile))
    assert ds.GetGeoTransform() == (727000, 4, 0, 6172000, 0, -4)
    assert ds.RasterXSize == 250
    assert ds.RasterYSize == 250
    band = ds.GetRasterBand(1)
    assert band.DataType == gdal.GDT_Float32
    assert band.GetNoDataValue() == -99
    ds = None

    outfile = tmp_path / "var.tif"
    ds = gdal.Open(str(outfile))
    assert ds.GetGeoTransform() == (727000, 4, 0, 6172000, 0, -4)
    assert ds.RasterXSize == 250
    assert ds.RasterYSize == 250
    band = ds.GetRasterBand(1)
    assert band.DataType == gdal.GDT_Float32
    assert band.GetNoDataValue() == -99
    ds = None 
Example 8
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions.py    MIT License 6 votes vote down vote up
def toRaster(array, path, geometry, srs, navalue=-9999):
    """
    Writes a single array to a raster with coordinate system and geometric
    information.

    path = target path
    srs = spatial reference system
    """
    xpixels = array.shape[1]
    ypixels = array.shape[0]
    path = path.encode('utf-8')
    image = gdal.GetDriverByName("GTiff").Create(path, xpixels, ypixels,
                                1, gdal.GDT_Float32)
    image.SetGeoTransform(geometry)
    image.SetProjection(srs)
    image.GetRasterBand(1).WriteArray(array)
    image.GetRasterBand(1).SetNoDataValue(navalue) 
Example 9
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions.py    MIT License 6 votes vote down vote up
def toRasters(arraylist, path, geometry, srs):
    """
    Writes a list of 2d arrays, or a 3d array, to a series of tif files.

    Arraylist format = [[name,array],[name,array],....]
    path = target path
    geometry = gdal geometry object
    srs = spatial reference system object
    """
    if path[-2:] == "\\":
        path = path
    else:
        path = path + "\\"
    sample = arraylist[0][1]
    ypixels = sample.shape[0]
    xpixels = sample.shape[1]
    for ray in  tqdm(arraylist):
        image = gdal.GetDriverByName("GTiff").Create(os.path.join(path,
                                                              ray[0] + ".tif"),
                                    xpixels, ypixels, 1, gdal.GDT_Float32)
        image.SetGeoTransform(geometry)
        image.SetProjection(srs)
        image.GetRasterBand(1).WriteArray(ray[1]) 
Example 10
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions_v3.py    MIT License 6 votes vote down vote up
def toRaster(array, path, geometry, srs, navalue=-9999):
    """
    Writes a single array to a raster with coordinate system and geometric
    information.

    path = target path
    srs = spatial reference system
    """
    xpixels = array.shape[1]    
    ypixels = array.shape[0]
    path = path.encode('utf-8')
    image = gdal.GetDriverByName("GTiff").Create(path, xpixels, ypixels,
                                1, gdal.GDT_Float32)
    image.SetGeoTransform(geometry)
    image.SetProjection(srs)
    image.GetRasterBand(1).WriteArray(array)
    image.GetRasterBand(1).SetNoDataValue(navalue) 
Example 11
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions_v3.py    MIT License 6 votes vote down vote up
def toRasters(arraylist, path, geometry, srs):
    """
    Writes a list of 2d arrays, or a 3d array, to a series of tif files.

    Arraylist format = [[name,array],[name,array],....]
    path = target path
    geometry = gdal geometry object
    srs = spatial reference system object
    """
    if path[-2:] == "\\":
        path = path
    else:
        path = path + "\\"
    sample = arraylist[0][1]
    ypixels = sample.shape[0]
    xpixels = sample.shape[1]
    for ray in  tqdm(arraylist):
        image = gdal.GetDriverByName("GTiff").Create(os.path.join(path,
                                                              ray[0] + ".tif"),
                                    xpixels, ypixels, 1, gdal.GDT_Float32)
        image.SetGeoTransform(geometry)
        image.SetProjection(srs)
        image.GetRasterBand(1).WriteArray(ray[1]) 
Example 12
Project: PyGeoC   Author: lreis2415   File: postTauDEM.py    MIT License 6 votes vote down vote up
def output_compressed_dinf(dinfflowang, compdinffile, weightfile):
        """Output compressed Dinf flow direction and weight to raster file
        Args:
            dinfflowang: Dinf flow direction raster file
            compdinffile: Compressed D8 flow code
            weightfile: The correspond weight
        """
        dinf_r = RasterUtilClass.read_raster(dinfflowang)
        data = dinf_r.data
        xsize = dinf_r.nCols
        ysize = dinf_r.nRows
        nodata_value = dinf_r.noDataValue

        cal_dir_code = frompyfunc(DinfUtil.compress_dinf, 2, 3)
        updated_angle, dir_code, weight = cal_dir_code(data, nodata_value)

        RasterUtilClass.write_gtiff_file(dinfflowang, ysize, xsize, updated_angle,
                                         dinf_r.geotrans, dinf_r.srs, DEFAULT_NODATA, GDT_Float32)
        RasterUtilClass.write_gtiff_file(compdinffile, ysize, xsize, dir_code,
                                         dinf_r.geotrans, dinf_r.srs, DEFAULT_NODATA, GDT_Int16)
        RasterUtilClass.write_gtiff_file(weightfile, ysize, xsize, weight,
                                         dinf_r.geotrans, dinf_r.srs, DEFAULT_NODATA, GDT_Float32) 
Example 13
Project: PyGeoC   Author: lreis2415   File: raster.py    MIT License 6 votes vote down vote up
def __init__(self, n_rows, n_cols, data, nodata_value=None, geotransform=None,
                 srs=None, datatype=GDT_Float32):
        """Constructor."""
        self.nRows = n_rows
        self.nCols = n_cols
        self.data = numpy.copy(data)
        self.noDataValue = nodata_value
        self.geotrans = geotransform
        self.srs = srs
        self.dataType = datatype

        self.dx = geotransform[1]
        self.xMin = geotransform[0]
        self.xMax = geotransform[0] + n_cols * geotransform[1]
        self.yMax = geotransform[3]
        self.yMin = geotransform[3] + n_rows * geotransform[5]
        self.validZone = self.data != self.noDataValue
        self.validValues = numpy.where(self.validZone, self.data, numpy.nan) 
Example 14
Project: PyGeoC   Author: lreis2415   File: raster.py    MIT License 6 votes vote down vote up
def raster_to_gtiff(tif, geotif, change_nodata=False, change_gdal_type=False):
        """Converting Raster format to GeoTIFF.

        Args:
            tif: source raster file path.
            geotif: output raster file path.
            change_nodata: change NoDataValue to -9999 or not.
            gdal_type (:obj:`pygeoc.raster.GDALDataType`): GDT_Float32 as default.
            change_gdal_type: If True, output the Float32 data type.
        """
        rst_file = RasterUtilClass.read_raster(tif)
        nodata = rst_file.noDataValue
        if change_nodata:
            if not MathClass.floatequal(rst_file.noDataValue, DEFAULT_NODATA):
                nodata = DEFAULT_NODATA
                nodata_array = numpy.ones((rst_file.nRows, rst_file.nCols)) * rst_file.noDataValue
                nodata_check = numpy.isclose(rst_file.data, nodata_array)
                rst_file.data[nodata_check] = DEFAULT_NODATA
                # rst_file.data[rst_file.data == rst_file.noDataValue] = DEFAULT_NODATA
        gdal_type = rst_file.dataType
        if change_gdal_type:
            gdal_type = GDT_Float32
        RasterUtilClass.write_gtiff_file(geotif, rst_file.nRows, rst_file.nCols, rst_file.data,
                                         rst_file.geotrans, rst_file.srs, nodata,
                                         gdal_type) 
Example 15
Project: pygeotools   Author: dshean   File: geolib.py    MIT License 6 votes vote down vote up
def mem_ds(res, extent, srs=None, dtype=gdal.GDT_Float32):
    """Create a new GDAL Dataset in memory

    Useful for various applications that require a Dataset
    """
    #These round down to int
    #dst_ns = int((extent[2] - extent[0])/res)
    #dst_nl = int((extent[3] - extent[1])/res)
    #This should pad by 1 pixel, but not if extent and res were calculated together to give whole int
    dst_ns = int((extent[2] - extent[0])/res + 0.99)
    dst_nl = int((extent[3] - extent[1])/res + 0.99)
    m_ds = gdal.GetDriverByName('MEM').Create('', dst_ns, dst_nl, 1, dtype)
    m_gt = [extent[0], res, 0, extent[3], 0, -res]
    m_ds.SetGeoTransform(m_gt)
    if srs is not None:
        m_ds.SetProjection(srs.ExportToWkt())
    return m_ds

#Modify proj/gt of dst_fn in place 
Example 16
Project: pygeotools   Author: dshean   File: geolib.py    MIT License 6 votes vote down vote up
def gdaldem_mem_ma(ma, ds=None, res=None, extent=None, srs=None, processing='hillshade', returnma=False, computeEdges=False):
    """
    Wrapper to allow gdaldem calculations for arbitrary NumPy masked array input
    Untested, work in progress placeholder
    Should only need to specify res, can caluclate local gt, cartesian srs
    """
    if ds is None:
        ds = mem_ds(res, extent, srs=None, dtype=gdal.GDT_Float32)
    else:
        ds = mem_ds_copy(ds)
    b = ds.GetRasterBand(1)
    b.WriteArray(ma)
    out = gdaldem_mem_ds(ds, processing=processing, returnma=returnma)
    return out 

#Should add gdal.DEMProcessingOptions support 
Example 17
Project: pygeotools   Author: dshean   File: malib.py    MIT License 6 votes vote down vote up
def write_stats(self):
        #if not hasattr(self, 'stack_count'):
        stat_list = ['stack_count', 'stack_mean', 'stack_std', 'stack_min', 'stack_max']
        if self.med:
            stat_list.append('stack_med')
            stat_list.append('stack_nmad')
        if any([not hasattr(self, i) for i in stat_list]):
            self.compute_stats()
        print("Writing out stats")
        #Create dummy ds - might want to use vrt here instead
        driver = gdal.GetDriverByName("MEM")
        ds = driver.Create('', self.ma_stack.shape[2], self.ma_stack.shape[1], 1, gdal.GDT_Float32)
        ds.SetGeoTransform(self.gt)
        ds.SetProjection(self.proj)
        #Write out with malib, should preserve ma type
        out_prefix = os.path.splitext(self.stack_fn)[0]
        iolib.writeGTiff(self.stack_count, out_prefix+'_count.tif', ds)
        iolib.writeGTiff(self.stack_mean, out_prefix+'_mean.tif', ds)
        iolib.writeGTiff(self.stack_std, out_prefix+'_std.tif', ds)
        iolib.writeGTiff(self.stack_min, out_prefix+'_min.tif', ds)
        iolib.writeGTiff(self.stack_max, out_prefix+'_max.tif', ds)
        if self.med:
            iolib.writeGTiff(self.stack_med, out_prefix+'_med.tif', ds)
            iolib.writeGTiff(self.stack_nmad, out_prefix+'_nmad.tif', ds) 
Example 18
Project: pygeotools   Author: dshean   File: malib.py    MIT License 6 votes vote down vote up
def write_trend(self):
        #stat_list = ['stack_trend', 'stack_intercept', 'stack_detrended_std', 'stack_rsquared']
        stat_list = ['stack_trend', 'stack_intercept', 'stack_detrended_std']
        if any([not hasattr(self, i) for i in stat_list]):
            self.compute_trend()
        if self.stack_trend.count() > 0:
            print("Writing out trend")
            #Create dummy ds - might want to use vrt here instead
            driver = gdal.GetDriverByName("MEM")
            ds = driver.Create('', self.ma_stack.shape[2], self.ma_stack.shape[1], 1, gdal.GDT_Float32)
            ds.SetGeoTransform(self.gt)
            ds.SetProjection(self.proj)
            #Write out with malib, should preserve ma type
            out_prefix = os.path.splitext(self.stack_fn)[0]
            iolib.writeGTiff(self.stack_trend, out_prefix+'_trend.tif', ds)
            iolib.writeGTiff(self.stack_intercept, out_prefix+'_intercept.tif', ds)
            iolib.writeGTiff(self.stack_detrended_std, out_prefix+'_detrended_std.tif', ds)
            #iolib.writeGTiff(self.stack_rsquared, out_prefix+'_rsquared.tif', ds) 
Example 19
Project: fire-data-processing   Author: ANU-WALD   File: update_flammability_mosaic.py    GNU Affero General Public License v3.0 6 votes vote down vote up
def update_flammability_mosaic(date, n_band, dst, tmp, comp):
            
    flam = compose_mosaic(date, n_band, "flammability", gdal.GDT_Float32)
    anom = compose_mosaic(date, n_band, "anomaly", gdal.GDT_Float32)
    qmask = compose_mosaic(date, n_band, "quality_mask", gdal.GDT_Byte)
    
    tmp_file = os.path.join(tmp, uuid.uuid4().hex + ".nc")
    d = datetime.utcfromtimestamp(date.astype('O')/1e9)
    pack_flammability_mosaic(d, flam, anom, qmask, tmp_file)

    if not os.path.isfile(dst):
        shutil.move(tmp_file, dst)
    else:
        tmp_file2 = os.path.join(tmp, uuid.uuid4().hex + ".nc")
        os.system("cdo mergetime {} {} {}".format(tmp_file, dst, tmp_file2))
        os.remove(tmp_file)
        if comp:
            # In case the UNLIMITED time dimension is causing problems. nccopy has: 
            # [-u] convert unlimited dimensions to fixed size in output
            tmp_file3 = os.path.join(tmp, uuid.uuid4().hex + ".nc")
            os.system("nccopy -d 4 -c 'time/4,y/240,x/240' {} {}".format(tmp_file2, tmp_file3))
            os.remove(tmp_file2)
            shutil.move(tmp_file3, dst)
        else: 
            shutil.move(tmp_file2, dst) 
Example 20
Project: fire-data-processing   Author: ANU-WALD   File: update_fmc_mosaic.py    GNU Affero General Public License v3.0 6 votes vote down vote up
def update_fmc_mosaic(date, n_band, dst, tmp, comp):
            
    fmc_mean = compose_mosaic(date, n_band, "lfmc_mean", gdal.GDT_Float32)
    fmc_stdv = compose_mosaic(date, n_band, "lfmc_stdv", gdal.GDT_Float32)
    q_mask = compose_mosaic(date, n_band, "quality_mask", gdal.GDT_Byte)
    
    tmp_file = os.path.join(tmp, uuid.uuid4().hex + ".nc")
    d = datetime.utcfromtimestamp(date.astype('O')/1e9)
    pack_fmc_mosaic(d, fmc_mean, fmc_stdv, q_mask, tmp_file)

    if not os.path.isfile(dst):
        shutil.move(tmp_file, dst)
    else:
        tmp_file2 = os.path.join(tmp, uuid.uuid4().hex + ".nc")
        os.system("cdo mergetime {} {} {}".format(tmp_file, dst, tmp_file2))
        os.remove(tmp_file)
        if comp:
            # In case the UNLIMITED time dimension is causing problems. nccopy has: 
            # [-u] convert unlimited dimensions to fixed size in output
            tmp_file3 = os.path.join(tmp, uuid.uuid4().hex + ".nc")
            os.system("nccopy -d 4 -c 'time/4,y/240,x/240' {} {}".format(tmp_file2, tmp_file3))
            os.remove(tmp_file2)
            shutil.move(tmp_file3, dst)
        else: 
            shutil.move(tmp_file2, dst) 
Example 21
Project: CrownSeg   Author: rachit-ranjan16   File: calculate_ndvi.py    GNU General Public License v3.0 6 votes vote down vote up
def createOutputImage(self, outFilename, inDataset):
        '''
            Creates Output Image
        '''
        # Define the image driver to be used
        # This defines the output file format (e.g., GeoTiff)
        driver = gdal.GetDriverByName("GTiff")
        # Check that this driver can create a new file.
        metadata = driver.GetMetadata()
        if gdal.DCAP_CREATE in metadata and metadata[gdal.DCAP_CREATE] != 'YES':
            log.error('Driver GTIFF does not support Create()')
            sys.exit(-1)
        # Get the spatial information from the input file
        geoTransform = inDataset.GetGeoTransform()
        geoProjection = inDataset.GetProjection()
        # Create an output file of the same size as the inputted
        # image but with only 1 output image band.
        newDataset = driver.Create(outFilename, inDataset.RasterXSize,
                                   inDataset.RasterYSize, 1, gdal.GDT_Float32)
        # Define the spatial information for the new image.
        newDataset.SetGeoTransform(geoTransform)
        newDataset.SetProjection(geoProjection)
        return newDataset 
Example 22
Project: pyBRAT   Author: Riverscapes   File: bdws.py    GNU General Public License v3.0 6 votes vote down vote up
def writeArrayToRaster(self, file, array, lowernd, uppernd):
        """
        Save numpy array as a GeoTiff raster concurrent with the input DEM.

        :param file: Path to save file.
        :param array: Numpy array of data.
        :param lowernd: Lowest data value.
        :param uppernd: Highest data value.

        :return: None
        """
        if array.shape == self.dem.shape:
            ds = self.driverTiff.Create(file, xsize=self.demDS.RasterXSize, ysize=self.demDS.RasterYSize, bands=1, eType=gdal.GDT_Float32)
            ds.SetGeoTransform(self.geot)
            ds.SetProjection(self.prj)
            array[np.isnan(array)] = -9999.0
            array[array < lowernd] = -9999.0
            array[array > uppernd] = -9999.0
            ds.GetRasterBand(1).WriteArray(array)
            ds.GetRasterBand(1).FlushCache()
            ds.GetRasterBand(1).SetNoDataValue(-9999.0)
            ds = None
        else:
            print "output shape different from input DEM" 
Example 23
Project: pyBRAT   Author: Riverscapes   File: bdflopy.py    GNU General Public License v3.0 6 votes vote down vote up
def createDatasets(self, filelist):
        """
        Create GDAL raster datasets.

        :param filelist: List of paths where raster datasets will be created.

        :return: List of GDAL datasets
        """
        datasetlist = []
        for file in filelist:
            #create raster dataset
            datasetlist.append(self.driver.Create(file, self.xsize, self.ysize, 1, gdal.GDT_Float32))
            #set projection to match input dem
            datasetlist[-1].SetProjection(self.prj)
            #set geotransform to match input dem
            datasetlist[-1].SetGeoTransform(self.geot)
        return datasetlist 
Example 24
Project: geophys_utils   Author: GeoscienceAustralia   File: _array2file.py    Apache License 2.0 5 votes vote down vote up
def array2file(data_arrays, 
               projection, 
               geotransform, 
               file_path, 
               file_format='GTiff', 
               dtype=gdal.GDT_Float32):
    '''
    Function to output array to any(?) GDAL supported raster format and return GDAL dataset
    Formats defined in http://www.gdal.org/formats_list.html
    '''
    data_array_shape = data_arrays[0].shape
    assert [data_array.shape for data_array in data_arrays].count(data_array_shape) == len(data_arrays), 'Data arrays are of different shape'

    driver = gdal.GetDriverByName(file_format)
    gdal_dataset = driver.Create(file_path, 
                                 data_array_shape[1], data_array_shape[0], # Array must be ordered yx
                                 len(data_arrays), 
                                 dtype)
    gdal_dataset.SetGeoTransform(geotransform)
    gdal_dataset.SetProjection(projection)
    
    # Write arrays to file
    for band_index in range(len(data_arrays)):
        raster_band = gdal_dataset.GetRasterBand(band_index+1)
        raster_band.WriteArray(data_arrays[band_index])
        raster_band.FlushCache()
        
    gdal_dataset.FlushCache()
    return gdal_dataset 
Example 25
Project: sentinelloader   Author: flaviostutz   File: utils.py    MIT License 5 votes vote down vote up
def saveGeoTiff(imageDataFloat, outputFile, geoTransform, projection):
    driver = gdal.GetDriverByName('GTiff')
    image_data = driver.Create(outputFile, imageDataFloat.shape[1], imageDataFloat.shape[0], 1, gdal.GDT_Float32)
    image_data.GetRasterBand(1).WriteArray(imageDataFloat)
    image_data.SetGeoTransform(geoTransform) 
    image_data.SetProjection(projection)
    image_data.FlushCache()
    image_data=None 
Example 26
Project: buzzard   Author: airware   File: tools.py    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 27
Project: buzzard   Author: airware   File: _gdal_gdt_conv.py    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 28
Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_GDALIO.py    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 29
Project: Landsat8LST_SWA   Author: eduard-kazakov   File: l8_lst_swa_common_lib.py    GNU General Public License v2.0 5 votes vote down vote up
def Landsat8_simple_temperature (mode, raster_path, channel_number, metadata_path, output_path=None, base_raster=None):
    # mode 0 - generate tiff
    # mode 1 - return np array

    """
    Return brigtness temperature for Landsat's 8 B10 or B11.
    :param mode: int; mode 0 - generate tiff; mode 1 - return np array
    :param raster_path: string; input raster path
    :param channel_number: int; number of channel to be processed (10 or 11)
    :param metadata_path: string; path to metadata file
    :param output_path: string; path for output brightness temperature file
    :param base_raster: base geotiff layer (needed if mode = 0)
    :return: string; path for output radiance file
    """
    try:
        int(channel_number)
    except:
        raise NameError('Channel number is not correct')

    if (channel_number < 10) or (channel_number > 11):
        raise NameError('Channel number is not correct')

    landsat_radiance_band_array = Landsat8_DN_to_radiance(1,raster_path,channel_number,metadata_path)
    K1_constant = float(read_metadata_parameter(metadata_path,'K1_CONSTANT_BAND_' + str(channel_number)))
    K2_constant = float(read_metadata_parameter(metadata_path,'K2_CONSTANT_BAND_' + str(channel_number)))

    landsat_temperature_array = (K2_constant / np.log((K1_constant/landsat_radiance_band_array)+1)) - 273.15
    del landsat_radiance_band_array
    if mode == 0:
        base_raster = gdal.Open(re.sub("^\s+|\n|\r|\s+$", '', raster_path))
        cols = base_raster.RasterXSize
        rows = base_raster.RasterYSize
        cell_type = gdal.GDT_Float32
        driver_name = 'GTiff'
        projection = base_raster.GetProjection()
        transform = base_raster.GetGeoTransform()
        save_nparray_as_raster(landsat_temperature_array,driver_name,cell_type,cols,rows,projection,transform,output_path)
        del landsat_temperature_array
    else:
        return landsat_temperature_array 
Example 30
Project: lidar   Author: giswqs   File: filling.py    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 31
Project: eo_tools   Author: neel9102   File: L2_LAC_SST_preprocessing.py    GNU General Public License v3.0 5 votes vote down vote up
def output_file(output_name,output_array):
  
    driver = gdal.GetDriverByName("GTiff")
    outDataset = driver.Create(output_name, ncols, nrows, 1, gdal.GDT_Float32)
    outDataset.SetGeoTransform(geo_transform_1)
    #srs = osr.SpatialReference()
    #srs.ImportFromEPSG(4326) #WGS84 lat long.
    outDataset.SetProjection(proj_1)
    
    outBand = outDataset.GetRasterBand(1)
    outBand.WriteArray(output_array,0,0)
    
    outBand.SetNoDataValue(np.nan)
    outBand.FlushCache() 
Example 32
Project: unmixing   Author: arthur-e   File: utils.py    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 33
Project: watools   Author: wateraccounting   File: data_conversions.py    Apache License 2.0 5 votes vote down vote up
def Save_as_tiff(name='', data='', geo='', projection=''):
    """
    This function save the array as a geotiff

    Keyword arguments:
    name -- string, directory name
    data -- [array], dataset of the geotiff
    geo -- [minimum lon, pixelsize, rotation, maximum lat, rotation,
            pixelsize], (geospatial dataset)
    projection -- integer, the EPSG code
    """
    # save as a geotiff
    driver = gdal.GetDriverByName("GTiff")
    dst_ds = driver.Create(name, int(data.shape[1]), int(data.shape[0]), 1,
                           gdal.GDT_Float32, ['COMPRESS=LZW'])
    srse = osr.SpatialReference()
    if projection == '':
        srse.SetWellKnownGeogCS("WGS84")

    else:
        try:
            if not srse.SetWellKnownGeogCS(projection) == 6:
                srse.SetWellKnownGeogCS(projection)
            else:
                try:
                    srse.ImportFromEPSG(int(projection))
                except:
                    srse.ImportFromWkt(projection)
        except:
            try:
                srse.ImportFromEPSG(int(projection))
            except:
                srse.ImportFromWkt(projection)

    dst_ds.SetProjection(srse.ExportToWkt())
    dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
    dst_ds.SetGeoTransform(geo)
    dst_ds.GetRasterBand(1).WriteArray(data)
    dst_ds = None
    return() 
Example 34
Project: pyIntertidalDEM   Author: jamal919   File: core.py    Apache License 2.0 5 votes vote down vote up
def to_geotiff(self, fname, dtype=gdal.GDT_Float32, epsg='auto'):
        '''
        Save band data to geotiff to location passed by `to` with datatype
        defined by `dtype`

        argument:
            fname: string
                The filename to be saved
            dtype: gdal data type
                Gdal datatype to be used for saving, default `gdal.GDT_Float32`
            epsg: epsg code
                epsg code to reproject the data. `auto` saves the data to
                original projection. Default `auto` (only option)

        '''
        row, col, nband = self.rgb.shape
        
        if epsg=='auto':
            driver = gdal.GetDriverByName('GTiff')
            gtiff = driver.Create(fname, row, col, nband, dtype)
            gtiff.GetRasterBand(1).WriteArray(self.rgb[:, :, 0])
            gtiff.GetRasterBand(2).WriteArray(self.rgb[:, :, 1])
            gtiff.GetRasterBand(3).WriteArray(self.rgb[:, :, 2])
            gtiff.SetGeoTransform(self.geotransform)
            gtiff.SetProjection(self.projection)
            gtiff.FlushCache()
            gtiff = None
        else:
            raise NotImplementedError 
Example 35
Project: Semantic_Segmentation_Keras   Author: liuph0119   File: predicting.py    Apache License 2.0 5 votes vote down vote up
def save_prediction(arr, dst_path, to_prob=False, with_geo=False, spatial_reference=None):
    """ save prediction.
    :param arr: 2D/3D array.
    :param dst_path: string.
    :param to_prob: bool, whether to save probability. If True, a *.npy or a float-formatted tiff file will be created.
    :param with_geo: bool, whether to save spatial reference.
    :param spatial_reference: string, or dict. If the data type is string, this must be a path of a existing tiff file,
        else, must be a dict contains "geotransform" and "projection"
    :return:
    """
    if with_geo:
        if type(spatial_reference) is str and os.path.exists(spatial_reference):
            img_info = get_image_info(spatial_reference, get_geotransform=True, get_projection=True)
        elif type(spatial_reference) is dict \
                and "geotransform" in spatial_reference and "projection" in spatial_reference:
            img_info = spatial_reference
        else:
            raise ValueError("Invalid 'spatial_reference': None. Expected to be a dict or a existing tiff image path.")
        if to_prob:
            save_to_image_gdal(arr, dst_path, datatype=gdal.GDT_Float32,
                               geoTransform=img_info["geotransform"], proj=img_info["projection"])
        else:
            save_to_image_gdal(arr, dst_path, datatype=gdal.GDT_Byte,
                               geoTransform=img_info["geotransform"], proj=img_info["projection"])
    else:
        if to_prob:
            if ".npy" not in dst_path:
                dst_path = dst_path + ".npy"
            np.save(dst_path, arr)
        else:
            save_to_image(arr, dst_path) 
Example 36
Project: PyGeoC   Author: lreis2415   File: raster.py    MIT License 5 votes vote down vote up
def raster_reclassify(srcfile, v_dict, dstfile, gdaltype=GDT_Float32):
        """Reclassify raster by given classifier dict.

        Args:
            srcfile: source raster file.
            v_dict: classifier dict.
            dstfile: destination file path.
            gdaltype (:obj:`pygeoc.raster.GDALDataType`): GDT_Float32 as default.
        """
        src_r = RasterUtilClass.read_raster(srcfile)
        src_data = src_r.data
        dst_data = numpy.copy(src_data)
        if gdaltype == GDT_Float32 and src_r.dataType != GDT_Float32:
            gdaltype = src_r.dataType
        no_data = src_r.noDataValue
        new_no_data = DEFAULT_NODATA
        if gdaltype in [GDT_Unknown, GDT_Byte, GDT_UInt16, GDT_UInt32]:
            new_no_data = 0
        if not MathClass.floatequal(new_no_data, src_r.noDataValue):
            if src_r.noDataValue not in v_dict:
                v_dict[src_r.noDataValue] = new_no_data
                no_data = new_no_data

        for (k, v) in iteritems(v_dict):
            dst_data[src_data == k] = v
        RasterUtilClass.write_gtiff_file(dstfile, src_r.nRows, src_r.nCols, dst_data,
                                         src_r.geotrans, src_r.srs, no_data, gdaltype) 
Example 37
Project: PyGeoC   Author: lreis2415   File: raster.py    MIT License 5 votes vote down vote up
def write_gtiff_file(f_name, n_rows, n_cols, data, geotransform, srs, nodata_value,
                         gdal_type=GDT_Float32):
        """Output Raster to GeoTiff format file.

        Args:
            f_name: output gtiff file name.
            n_rows: Row count.
            n_cols: Col count.
            data: 2D array data.
            geotransform: geographic transformation.
            srs: coordinate system.
            nodata_value: nodata value.
            gdal_type (:obj:`pygeoc.raster.GDALDataType`): output raster data type,
                                                                  GDT_Float32 as default.
        """
        UtilClass.mkdir(os.path.dirname(FileClass.get_file_fullpath(f_name)))
        driver = gdal_GetDriverByName(str('GTiff'))
        try:
            ds = driver.Create(f_name, n_cols, n_rows, 1, gdal_type)
        except Exception:
            print('Cannot create output file %s' % f_name)
            return
        ds.SetGeoTransform(geotransform)
        try:
            ds.SetProjection(srs.ExportToWkt())
        except AttributeError or Exception:
            ds.SetProjection(srs)
        ds.GetRasterBand(1).SetNoDataValue(nodata_value)
        # if data contains numpy.nan, then replaced by nodata_value
        if isinstance(data, numpy.ndarray) and data.dtype in [numpy.dtype('int'),
                                                              numpy.dtype('float')]:
            data = numpy.where(numpy.isnan(data), nodata_value, data)
        ds.GetRasterBand(1).WriteArray(data)
        ds = None 
Example 38
Project: geobricks_qgis_plugin_trmm   Author: geobricks   File: conversions.py    GNU General Public License v2.0 5 votes vote down vote up
def Float32(dataset_or_band):
    return ConvertedDataset(dataset_or_band, gdal.GDT_Float32) 
Example 39
Project: dsarsim   Author: adamoferro   File: image.py    GNU General Public License v3.0 5 votes vote down vote up
def __init__(self,size=(0,0),ps=(0,0),value=None,nodatav=-9999):
        
        gdal.UseExceptions()
        
        self.matrix_type=numpy.float32
        self.file_type=gdal.GDT_Float32

        self.image=None
        self.size=size
        self.pixel_size=ps
        self.nodatav=nodatav

        if value is not None:
            self.empty(value) 
Example 40
Project: pygeotools   Author: dshean   File: malib.py    MIT License 5 votes vote down vote up
def write_datestack(self):
        #stat_list = ['dt_stack_ptp', 'dt_stack_mean', 'dt_stack_min', 'dt_stack_max', 'dt_stack_center']
        stat_list = ['dt_stack_ptp', 'dt_stack_min', 'dt_stack_max', 'dt_stack_center']
        if any([not hasattr(self, i) for i in stat_list]):
            #self.make_datestack()
            self.compute_dt_stats()
        print("Writing out datestack stats")
        #Create dummy ds - might want to use vrt here instead
        driver = gdal.GetDriverByName("MEM")
        ds = driver.Create('', self.dt_stack_ptp.shape[1], self.dt_stack_ptp.shape[0], 1, gdal.GDT_Float32)
        ds.SetGeoTransform(self.gt)
        ds.SetProjection(self.proj)
        #Write out with malib, should preserve ma type
        out_prefix = os.path.splitext(self.stack_fn)[0]
        iolib.writeGTiff(self.dt_stack_ptp, out_prefix+'_dt_ptp.tif', ds)
        self.dt_stack_ptp.set_fill_value(-9999)
        #iolib.writeGTiff(self.dt_stack_mean, out_prefix+'_dt_mean.tif', ds)
        #self.dt_stack_mean.set_fill_value(-9999)
        iolib.writeGTiff(self.dt_stack_min, out_prefix+'_dt_min.tif', ds)
        self.dt_stack_min.set_fill_value(-9999)
        iolib.writeGTiff(self.dt_stack_max, out_prefix+'_dt_max.tif', ds)
        self.dt_stack_max.set_fill_value(-9999)
        iolib.writeGTiff(self.dt_stack_center, out_prefix+'_dt_center.tif', ds)
        self.dt_stack_center.set_fill_value(-9999)

    #Note: want ot change the variable names min/max here
    #Might be better to save out as multiband GTiff here 
Example 41
Project: firedpy   Author: earthlab   File: functions.py    MIT License 5 votes vote down vote up
def rasterize(src, dst, attribute, resolution, crs, extent, all_touch=False,
              na=-9999):

    # Open shapefile, retrieve the layer
    src_data = ogr.Open(src)
    layer = src_data.GetLayer()

    # Use transform to derive coordinates and dimensions
    xmin = extent[0]
    ymin = extent[1]
    xmax = extent[2]
    ymax = extent[3]

    # Create the target raster layer
    cols = int((xmax - xmin)/resolution)
    rows = int((ymax - ymin)/resolution) + 1
    trgt = gdal.GetDriverByName("GTiff").Create(dst, cols, rows, 1,
                                gdal.GDT_Float32)
    trgt.SetGeoTransform((xmin, resolution, 0, ymax, 0, -resolution))

    # Add crs
    refs = osr.SpatialReference()
    refs.ImportFromWkt(crs)
    trgt.SetProjection(refs.ExportToWkt())

    # Set no value
    band = trgt.GetRasterBand(1)
    band.SetNoDataValue(na)

    # Set options
    if all_touch is True:
        ops = ["-at", "ATTRIBUTE=" + attribute]
    else:
        ops = ["ATTRIBUTE=" + attribute]

    # Finally rasterize
    gdal.RasterizeLayer(trgt, [1], layer, options=ops)

    # Close target an source rasters
    del trgt
    del src_data 
Example 42
Project: fire-data-processing   Author: ANU-WALD   File: update_flammability_mosaic.py    GNU Affero General Public License v3.0 5 votes vote down vote up
def compose_mosaic(date, n_band, var_name, data_type):
    d = datetime.utcfromtimestamp(date.astype('O')/1e9)
    lat0 = -10.
    lat1 = -44.
    lon0 = 113.
    lon1 = 154.
    res = 0.005

    x_size = int((lon1 - lon0)/res)
    y_size = int((lat1 - lat0)/(-1*res))
    lats = np.linspace(lat0, lat1+res, num=y_size)
    lons = np.linspace(lon0, lon1-res, num=x_size)
    
    src = gdal.GetDriverByName('MEM').Create('', tile_size, tile_size, 1, data_type,)

    geot = [lon0, res, 0., lat0, 0., -1*res]
    dst = gdal.GetDriverByName('MEM').Create('', x_size, y_size, 1, data_type,)
    if data_type == gdal.GDT_Float32:
        dst.GetRasterBand(1).WriteArray(np.ones((y_size, x_size), dtype=np.float32) * -9999.9)
    dst.SetGeoTransform(geot)
    dst.SetProjection(wgs84_wkt)

    for au_tile in au_tiles:
        fname = flam_stack_path.format(d.year, au_tile)
        stack = gdal.Open('NETCDF:"{}":{}'.format(fname, var_name))
        if stack is None:
            continue
        
        band = stack.GetRasterBand(n_band).ReadAsArray()
        src.GetRasterBand(1).WriteArray(band)
    
        src.SetGeoTransform(stack.GetGeoTransform())
        src.SetProjection(stack.GetProjection())
        err = gdal.ReprojectImage(src, dst, None, None, gdal.GRA_NearestNeighbour)
    
    return dst.ReadAsArray() 
Example 43
Project: fire-data-processing   Author: ANU-WALD   File: update_fmc_mosaic.py    GNU Affero General Public License v3.0 5 votes vote down vote up
def compose_mosaic(date, n_band, var_name, data_type):
    d = datetime.utcfromtimestamp(date.astype('O')/1e9)
    lat0 = -10.
    lat1 = -44.
    lon0 = 113.
    lon1 = 154.
    res = 0.005

    x_size = int((lon1 - lon0)/res)
    y_size = int((lat1 - lat0)/(-1*res))
    lats = np.linspace(lat0, lat1+res, num=y_size)
    lons = np.linspace(lon0, lon1-res, num=x_size)
    
    src = gdal.GetDriverByName('MEM').Create('', tile_size, tile_size, 1, data_type,)

    geot = [lon0, res, 0., lat0, 0., -1*res]
    dst = gdal.GetDriverByName('MEM').Create('', x_size, y_size, 1, data_type,)
    if data_type == gdal.GDT_Float32:
        dst.GetRasterBand(1).WriteArray(np.ones((y_size, x_size), dtype=np.float32) * -9999.9)
    dst.SetGeoTransform(geot)
    dst.SetProjection(wgs84_wkt)

    for au_tile in au_tiles:
        fname = fmc_stack_path.format(d.year, au_tile)
        stack = gdal.Open('NETCDF:"{}":{}'.format(fname, var_name))
        if stack is None:
            continue
        
        band = stack.GetRasterBand(n_band).ReadAsArray()
        src.GetRasterBand(1).WriteArray(band)
    
        src.SetGeoTransform(stack.GetGeoTransform())
        src.SetProjection(stack.GetProjection())
        err = gdal.ReprojectImage(src, dst, None, None, gdal.GRA_NearestNeighbour)
    
    return dst.ReadAsArray() 
Example 44
Project: hrds   Author: EnvModellingGroup   File: raster_buffer.py    GNU General Public License v3.0 5 votes vote down vote up
def __write_raster__(self, filename, array, dx, origin, proj):
        """
        Write raster to file

        Args:
            filename: filename to write to
            array: the data. 2D numpy array
            dx: resolution (length 2 list)
            origin: the LLC coordinates
            proj: Projection space for the raster (wkt)

        Returns:
            Nothing
        """
        dst_filename = filename
        x_pixels = array.shape[1]
        y_pixels = array.shape[0]
        x_min = origin[0][0]
        y_max = origin[0][1] + dx[1]*y_pixels
        wkt_projection = proj

        driver = gdal.GetDriverByName('GTiff')
        dataset = driver.Create(
            dst_filename,
            x_pixels,
            y_pixels,
            1,
            gdal.GDT_Float32)

        dataset.SetGeoTransform((
            x_min,       # 0
            dx[1],       # 1
            0,           # 2
            y_max,       # 3
            0,           # 4
            -dx[0]))

        dataset.SetProjection(wkt_projection)
        dataset.GetRasterBand(1).WriteArray(array)
        dataset.FlushCache()
        return 
Example 45
Project: Hydro-DS   Author: CI-WATER   File: netcdfFunctions.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def project_and_resample_Array(input_array, srs_geotrs, srs_proj, Nxin, Nyin, reference_netcdf):  #, output_array):

    #srs_data = gdal.Open(input_raster, GA_ReadOnly)
    #srs_proj = srs_data.GetProjection() #osr.SpatialReference(wkt
    srs_data = gdal.GetDriverByName('MEM').Create('', Nxin, Nyin, 1, gdal.GDT_Float32)
    srs_data.SetGeoTransform(srs_geotrs)
    srs_data.SetProjection(srs_proj)
    srsband = srs_data.GetRasterBand(1)
    srsband.WriteArray(input_array)
    srsband.FlushCache()

    ref_data = gdal.Open(reference_netcdf, GA_ReadOnly)
    ref_proj = ref_data.GetProjection()
    ref_geotrs = ref_data.GetGeoTransform()
    Ncols = ref_data.RasterXSize
    Nrows = ref_data.RasterYSize
    ref_data = None

    out_data = gdal.GetDriverByName('MEM').Create('', Ncols, Nrows, 1, gdal.GDT_Float32)
    out_data.SetGeoTransform(ref_geotrs)
    out_data.SetProjection(ref_proj)

    gdal.ReprojectImage(srs_data,out_data,srs_proj,ref_proj, gdal.GRA_Bilinear )
    output_array = out_data.ReadAsArray()

    srs_data = None
    out_data = None
    return output_array 
Example 46
Project: pyBRAT   Author: Riverscapes   File: bdws.py    GNU General Public License v3.0 5 votes vote down vote up
def writeDamLocationRaster(self):
        """
        Write dam locations to raster.

        :return: None
        """
        ds = self.driverTiff.Create(self.outDir + "/damID.tif", xsize=self.demDS.RasterXSize, ysize=self.demDS.RasterYSize, bands=1,
                                    eType=gdal.GDT_Float32)
        ds.SetGeoTransform(self.geot)
        ds.SetProjection(self.prj)
        ds.GetRasterBand(1).WriteArray(self.idOut)
        ds.GetRasterBand(1).FlushCache()
        ds.GetRasterBand(1).SetNoDataValue(-9999.0)
        ds = None 
Example 47
Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_GDALIO.py    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 48
Project: Landsat8LST_SWA   Author: eduard-kazakov   File: l8_lst_swa_common_lib.py    GNU General Public License v2.0 4 votes vote down vote up
def Landsat8_DN_to_radiance (mode, raster_path, channel_number, metadata_path, output_path = None):
    """
    Convert raw Landsat8 scene from DN to radiance
    :param mode: int; mode 0 - generate tiff; mode 1 - return np array
    :param raster_path: string; input raster path
    :param channel_number: int; number of channel to be processed
    :param metadata_path: string; path to metadata file
    :param output_path: string; path for output radiance file
    :return:
    """
    try:
        int(channel_number)
    except:
        raise NameError('Channel number is not correct')

    if (int(channel_number)) < 1 or (int(channel_number) > 11):
        raise NameError('Channel number is not correct')

    metadata_channel_rad_maximum_str = 'RADIANCE_MAXIMUM_BAND_' + str(channel_number)
    metadata_channel_rad_minimum_str = 'RADIANCE_MINIMUM_BAND_' + str(channel_number)

    metadata_channel_quantize_max_str = 'QUANTIZE_CAL_MAX_BAND_' + str(channel_number)
    metadata_channel_quantize_min_str = 'QUANTIZE_CAL_MIN_BAND_' + str(channel_number)

    rad_maximum = float(read_metadata_parameter(metadata_path, metadata_channel_rad_maximum_str))
    rad_minimum = float(read_metadata_parameter(metadata_path, metadata_channel_rad_minimum_str))
    quantize_maximum = float(read_metadata_parameter(metadata_path, metadata_channel_quantize_max_str))
    quantize_minimum = float(read_metadata_parameter(metadata_path, metadata_channel_quantize_min_str))

    #print raster_path.replace('/','\\')
    landsat_dn_band = gdal.Open(re.sub("^\s+|\n|\r|\s+$", '', raster_path))
    #print landsat_dn_band
    landsat_dn_band_array = np.array(landsat_dn_band.GetRasterBand(1).ReadAsArray().astype(np.float32))
    print 'a???'
    landsat_radiance_band_array = ((rad_maximum - rad_minimum)/(quantize_maximum - quantize_minimum))*(landsat_dn_band_array - quantize_minimum) + rad_minimum
    print 'b???'
    if mode == 0:
        # Write result
        cols = landsat_dn_band.RasterXSize
        rows = landsat_dn_band.RasterYSize
        cell_type = gdal.GDT_Float32
        driver_name = 'GTiff'
        projection = landsat_dn_band.GetProjection()
        transform = landsat_dn_band.GetGeoTransform()
        save_nparray_as_raster(landsat_radiance_band_array,driver_name,cell_type,cols,rows,projection,transform,output_path)
        del landsat_radiance_band_array
        return
    else:
        return landsat_radiance_band_array 
Example 49
Project: Landsat8LST_SWA   Author: eduard-kazakov   File: l8_lst_swa_main.py    GNU General Public License v2.0 4 votes vote down vote up
def getInputs(self):
        tempPath = os.path.dirname(os.path.abspath(__file__)) + '\\temp\\'
        if self.ui.satTabDataTypeRawRadioButton.isChecked():
            metadataPath = self.ui.satTabMTLPathLine.text()
            metadataBasePath = os.path.dirname(metadataPath)
            mtl_dict, B10BrightnessTemperature, B11BrightnessTemperature, LSE10Array, LSE11Array = self.prepareRawLandsat()
            #print 'lol'

        if self.ui.waterVaporGRIDRadioButton.isChecked():
            #print 'adj'
            baseRasterPath = re.sub("^\s+|\n|\r|\s+$", '', metadataBasePath + '/' + mtl_dict['mtl_band10'])
            baseRaster = QgsRasterLayer(baseRasterPath,'B10 Base Layer')
            #print baseRaster.isValid()
            #print baseRaster.extent()
            waterVaporSourceRaster = l8_lst_swa_common_lib.getLayerByName(self.ui.waterVaporGRIDComboBox.currentText())
            #print waterVaporSourceRaster.isValid()
            l8_lst_swa_common_lib.adjustRasterToBaseRaster(baseRaster,waterVaporSourceRaster,tempPath + 'waterVapor.tif')
            waterVapor = gdal.Open (tempPath + 'waterVapor.tif')
            waterVaporArray = np.array(waterVapor.GetRasterBand(1).ReadAsArray().astype(np.float32))

        print waterVaporArray.shape
        print B10BrightnessTemperature.shape
        print B11BrightnessTemperature.shape
        print LSE10Array.shape
        print LSE11Array.shape

        #LSTArray = l8_lst_swa_core.getLSTWithSWAForArray(waterVaporArray,LSE10Array,LSE11Array,B10BrightnessTemperature,B11BrightnessTemperature)
        LSTArray = (B10BrightnessTemperature + B11BrightnessTemperature) / 2
        cols = waterVapor.RasterXSize
        rows = waterVapor.RasterYSize
        cell_type = gdal.GDT_Float32
        driver_name = 'GTiff'
        projection = waterVapor.GetProjection()
        transform = waterVapor.GetGeoTransform()
#
        output_path = self.ui.outputLSTLine.text()
#
        l8_lst_swa_common_lib.save_nparray_as_raster(LSTArray,driver_name,cell_type,cols,rows,projection,transform,output_path)

        del waterVaporArray
        del B10BrightnessTemperature
        del B11BrightnessTemperature
        del LSE10Array
        del LSE11Array
        gc.collect() 
Example 50
Project: watools   Author: wateraccounting   File: raster_conversions.py    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)