Python osgeo.gdal.GA_Update() Examples

The following are code examples for showing how to use osgeo.gdal.GA_Update(). 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: arsf_dem_scripts   Author: pmlrsg   File: dem_utilities.py    GNU General Public License v3.0 6 votes vote down vote up
def set_nodata_value(in_file, nodata_value):
    """
    Sets nodata value for the first band of a GDAL dataset.

    Arguments:

    * in_file - path to existing GDAL dataset.
    * nodata_value - nodata value

    """

    if not HAVE_GDAL:
        raise ImportError('Could not import GDAL')

    gdal_ds = gdal.Open(in_file, gdal.GA_Update)
    gdal_ds.GetRasterBand(1).SetNoDataValue(nodata_value)
    gdal_ds = None 
Example 2
Project: pygeotools   Author: dshean   File: geolib.py    MIT License 6 votes vote down vote up
def copyproj(src_fn, dst_fn, gt=True):
    """Copy projection and geotransform from one raster file to another
    """
    src_ds = gdal.Open(src_fn, gdal.GA_ReadOnly)
    dst_ds = gdal.Open(dst_fn, gdal.GA_Update)
    dst_ds.SetProjection(src_ds.GetProjection())
    if gt:
        src_gt = np.array(src_ds.GetGeoTransform())
        src_dim = np.array([src_ds.RasterXSize, src_ds.RasterYSize])
        dst_dim = np.array([dst_ds.RasterXSize, dst_ds.RasterYSize])
        #This preserves dst_fn resolution
        if np.any(src_dim != dst_dim):
            res_factor = src_dim/dst_dim.astype(float)
            src_gt[[1, 5]] *= max(res_factor)
            #src_gt[[1, 5]] *= min(res_factor)
            #src_gt[[1, 5]] *= res_factor
        dst_ds.SetGeoTransform(src_gt)
    src_ds = None
    dst_ds = None 
Example 3
Project: soilscape_upscaling   Author: MiXIL   File: stack_bands.py    GNU General Public License v3.0 6 votes vote down vote up
def set_band_names(input_image, band_names_list, print_names=False):
    """
    A utility function to set band names.
    Taken from RSGISLib (rsgislib.org)

    Requires:

    * input_image
    * band_names_list - list of names for each band

    """
    dataset = gdal.Open(input_image, gdal.GA_Update)

    for i, band_name in enumerate(band_names_list):
        band = i+1
        img_band = dataset.GetRasterBand(band)
        # Check the image band is available
        if img_band is not None:
            if print_names:
                print('Setting Band {0} to "{1}"'.format(band, band_name))
            img_band.SetDescription(band_name)
        else:
            raise Exception("Could not open the image band: {}".format(band))

    dataset = None 
Example 4
Project: AtmosphericCorrection   Author: y-iikura   File: apm_util.py    MIT License 5 votes vote down vote up
def read_tif(fname):
  src = gdal.Open(fname, gdal.GA_Update)
  pdem = src.GetRasterBand(1)
  gt = src.GetGeoTransform()
  image = pdem.ReadAsArray()
  return [gt,image] 
Example 5
Project: pygeotools   Author: dshean   File: iolib.py    MIT License 5 votes vote down vote up
def set_ndv(dst_fn, ndv):
    dst_ds = gdal.Open(dst_fn, gdal.GA_Update)
    for n in range(1, dst_ds.RasterCount+1):
        b = dst_ds.GetRasterBand(1)
        b.SetNoDataValue(ndv)
    dst_ds = None

#Should overload these functions to handle fn, ds, or b
#Perhaps abstract, as many functions will need this functionality 
Example 6
Project: HistoricalMap   Author: nkarasiak   File: function_historical_map.py    GNU General Public License v2.0 5 votes vote down vote up
def reclassAndFillHole(self,rasterTemp,inClassNumber):
        """ [email protected] Reclass and file hole (with ndimage morphology binary fill hole) for last treatment (need raster and class number)"""
        dst_ds = gdal.Open(rasterTemp,gdal.GA_Update)
        srcband = dst_ds.GetRasterBand(1).ReadAsArray()
        # All data which is not forest is set to 0, so we fill all for the forest only, because it's a binary fill holes.            
        # Set selected class as 1                   
        srcband[srcband != inClassNumber]=0
        srcband[srcband == inClassNumber]=1
        
        srcband = ndimage.morphology.binary_fill_holes(srcband)
        
        dst_ds.GetRasterBand(1).WriteArray(srcband)
        dst_ds.FlushCache()
        return rasterTemp 
Example 7
Project: CoastSat   Author: kvos   File: SDS_tools.py    GNU General Public License v3.0 5 votes vote down vote up
def mask_raster(fn, mask):
    """
    Masks a .tif raster using GDAL.
    
    Arguments:
    -----------
        fn: str
            filepath + filename of the .tif raster
        mask: np.array
            array of boolean where True indicates the pixels that are to be masked
        
    Returns:    
    -----------
    overwrites the .tif file directly
        
    """ 
    
    # open raster
    raster = gdal.Open(fn, gdal.GA_Update)
    # mask raster
    for i in range(raster.RasterCount):
        out_band = raster.GetRasterBand(i+1)
        out_data = out_band.ReadAsArray()
        out_band.SetNoDataValue(0)
        no_data_value = out_band.GetNoDataValue()
        out_data[mask] = no_data_value
        out_band.WriteArray(out_data)
    # close dataset and flush cache
    raster = None


###################################################################################################
# UTILITIES
################################################################################################### 
Example 8
Project: arsf_dem_scripts   Author: pmlrsg   File: dem_utilities.py    GNU General Public License v3.0 4 votes vote down vote up
def add_dem_metadata(dem_name, dem_source=None, dem_filename=None,
                               other_items=None):
    """
    Adds metadata to DEM with original DEM source

    Arguments:

    * dem_name - Name of DEM file
    * dem_source (e.g., ASTER)
    * dem_filename (name of mosaic DEM was derived from)
    """
    if not HAVE_GDAL:
        raise ImportError('Could not import GDAL')

    # If ENVI dataset, write key to header file
    if get_gdal_type_from_path(dem_name) == 'ENVI':
        # Get ENVI header (using GDAL filelist function)
        dem_dataset  = gdal.Open(dem_name, gdal.GA_ReadOnly)
        if dem_dataset is None:
            raise Exception('Could not open {}'.format(dem_name))
        dem_header = dem_dataset.GetFileList()[1]
        dem_dataset = None

        # Open header and append text to bottom
        with open(dem_header,'a') as f:
            if dem_source is not None:
                f.write(';DEM Source={}\n'.format(dem_source))
            if dem_filename is not None:
                f.write(';Original DEM filename={}\n'.format(dem_filename))
            if other_items is not None:
                if type(other_items) is not dict:
                    raise TypeError('"other_items" must be a dictionary')
                for key in other_items.keys():
                    f.write(';{}={}\n'.format(key, other_items[key]))

    # If not try to write using GDAL metadata function.
    else:
        dem_dataset  = gdal.Open(dem_name, gdal.GA_Update)
        if dem_dataset is None:
            raise Exception('Could not open {} using GDAL to write metadata to'.format(dem_name))

        metadata = dem_dataset.GetMetadata()

        if dem_source is not None:
            metadata['DEM Source'] = dem_source
        if dem_filename is not None:
            metadata['Original DEM filename'] = dem_filename
        if other_items is not None:
            if type(other_items) is not dict:
                raise TypeError('"other_items" must be a dictionary')
            metadata.update(other_items)

        # Save updated metadata dictionary
        dem_dataset.SetMetadata(metadata)

        # Close dataset
        dem_dataset = None 
Example 9
Project: dask-geomodeling   Author: nens   File: utils.py    BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
def create_dataset(array, geo_transform=None, projection=None, no_data_value=None):
    """
    Create and return a gdal dataset.

    :param array: A numpy array.
    :param geo_transform: 6-tuple of floats
    :param projection: wkt projection string
    :param no_data_value: integer or float

    This is the fastest way to get a gdal dataset from a numpy array, but
    keep a reference to the array around, or a segfault will occur. Also,
    don't forget to call FlushCache() on the dataset after any operation
    that affects the array.
    """
    # prepare dataset name pointing to array
    datapointer = array.ctypes.data
    bands, lines, pixels = array.shape
    datatypecode = gdal_array.NumericTypeCodeToGDALTypeCode(array.dtype.type)
    datatype = gdal.GetDataTypeName(datatypecode)
    bandoffset, lineoffset, pixeloffset = array.strides
    # if projection is wrong there will be a segfault in RasterizeLayer
    projection = osr.GetUserInputAsWKT(str(projection))

    dataset_name_template = (
        "MEM:::"
        "DATAPOINTER={datapointer},"
        "PIXELS={pixels},"
        "LINES={lines},"
        "BANDS={bands},"
        "DATATYPE={datatype},"
        "PIXELOFFSET={pixeloffset},"
        "LINEOFFSET={lineoffset},"
        "BANDOFFSET={bandoffset}"
    )
    dataset_name = dataset_name_template.format(
        datapointer=datapointer,
        pixels=pixels,
        lines=lines,
        bands=bands,
        datatype=datatype,
        pixeloffset=pixeloffset,
        lineoffset=lineoffset,
        bandoffset=bandoffset,
    )

    # access the array memory as gdal dataset
    dataset = gdal.Open(dataset_name, gdal.GA_Update)

    # set additional properties from kwargs
    if geo_transform is not None:
        dataset.SetGeoTransform(geo_transform)
    if projection is not None:
        dataset.SetProjection(projection)
    if no_data_value is not None:
        for i in range(len(array)):
            dataset.GetRasterBand(i + 1).SetNoDataValue(no_data_value)

    return dataset