Python gdal.Open() Examples

def save_raster_simple(array, path, dst_filename):

    """ Save an array base on an existing raster """

    example = gdal.Open(path)
    x_pixels = array.shape[1]  # number of pixels in x
    y_pixels = array.shape[0]  # number of pixels in y
    bands = 1
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(dst_filename,x_pixels, y_pixels, bands, 
                            gdal.GDT_Int32)

    geotrans=example.GetGeoTransform()  #get GeoTranform from existed 'data0'
    proj=example.GetProjection() #you can get from a exsited tif or import 
    dataset.SetGeoTransform(geotrans)
    dataset.SetProjection(proj)

    dataset.GetRasterBand(1).WriteArray(array[:,:])

    dataset.FlushCache() 

def clip_raster_to_intersection(raster_to_clip_path, extent_raster_path, out_raster_path, is_landsat=False):
    """
    Clips one raster to the extent proivded by the other raster, and saves the result at out_raster_path.
    Assumes both raster_to_clip and extent_raster are in the same projection.
    Parameters
    ----------
    raster_to_clip_path
        The location of the raster to be clipped.
    extent_raster_path
        The location of the raster that will provide the extent to clip to
    out_raster_path
        A location for the finished raster
    """

    with TemporaryDirectory() as td:
        temp_aoi_path = os.path.join(td, "temp_clip.shp")
        get_extent_as_shp(extent_raster_path, temp_aoi_path)
        ext_ras = gdal.Open(extent_raster_path)
        proj = osr.SpatialReference(wkt=ext_ras.GetProjection())
        srs_id = int(proj.GetAttrValue('AUTHORITY', 1))
        clip_raster(raster_to_clip_path, temp_aoi_path, out_raster_path, srs_id, flip_x_y = is_landsat) 

def save_raster_memory(array, path):

    """ Save a raster into memory """

    example = gdal.Open(path)
    x_pixels = array.shape[1]  # number of pixels in x
    y_pixels = array.shape[0]  # number of pixels in y
    driver = gdal.GetDriverByName('MEM')
    dataset = driver.Create('',x_pixels, y_pixels, 1,gdal.GDT_Int32)
    dataset.GetRasterBand(1).WriteArray(array[:,:])

    # follow code is adding GeoTranform and Projection
    geotrans=example.GetGeoTransform()  #get GeoTranform from existed 'data0'
    proj=example.GetProjection() #you can get from a exsited tif or import 
    dataset.SetGeoTransform(geotrans)
    dataset.SetProjection(proj)

    return dataset 

def save_raster(array, path, dst_filename):
    """ Save the final multiband array based on an existing raster """

    example = gdal.Open(path)
    x_pixels = array.shape[2]  # number of pixels in x
    y_pixels = array.shape[1]  # number of pixels in y
    bands = array.shape[0]
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(dst_filename,x_pixels, 
                            y_pixels, bands ,gdal.GDT_Int32)

    geotrans=example.GetGeoTransform()  #get GeoTranform from existed 'data0'
    proj=example.GetProjection() #you can get from a exsited tif or import 
    dataset.SetGeoTransform(geotrans)
    dataset.SetProjection(proj)

    for b in range(bands):
        dataset.GetRasterBand(b+1).WriteArray(array[b,:,:])

    dataset.FlushCache() 

def test_mask_combination():
    os.chdir(os.path.dirname(os.path.abspath(__file__)))
    masks = [r"test_data/"+mask for mask in os.listdir("test_data") if mask.endswith(".msk")]
    try:
        os.remove("test_outputs/union_or_combination.tif")
        os.remove("test_outputs/intersection_and_combination.tif")
    except FileNotFoundError:
        pass
    pyeo.raster_manipulation.combine_masks(masks, "test_outputs/union_or_combination.tif",
                                           geometry_func="union", combination_func="or")
    pyeo.raster_manipulation.combine_masks(masks, "test_outputs/intersection_and_combination.tif",
                                           geometry_func="intersect", combination_func="and")
    mask_1 = gdal.Open("test_outputs/union_or_combination.tif")
    assert not mask_1.GetVirtualMemArray().all == False
    mask_2 = gdal.Open("test_outputs/intersection_and_combination.tif")
    assert not mask_2.GetVirtualMemArray().all == False 

def reproject_modis_to_geotiff(input_folder, dest_prj=32643):
    '''
    Function to convert all modis hdf in folder to geotiff using gdal
    required libraries: osr, gdal
    :param input_folder: where hds files are stored
    :param dest_prj: default is wgs 84 / UTM 43N (EPSG 32643)
    '''
    files_list = os.listdir(input_folder)
    for item in files_list:
        if fnmatch.fnmatch(item, '*.hdf'):
            input_file_name = input_folder + '/' + item
            output_file_name = input_folder + '/' + item[0:23]
            img = gdal.Open(input_file_name)
            subset = img.GetSubDatasets()
            in_raster = subset[0][0]
            reproject_dataset(dataset=in_raster, output_file_name=output_file_name+ '.tif')

# reproject
# reproject_modis_to_geotiff(input_folder="/media/kiruba/New Volume/MODIS/ET/scratch")

# This function will convert the rasterized clipper shapefile
# to a mask for use within GDAL. 

def fix_alpha_channel(filename):
    """
    Imagemagick's 'convert' command converts all channels to 8-bits, including the alpha channel.
    Unfortunately it converts the alpha value 255 to 1 when it does this. This is a hack to restore
    full transparency values to 255 in the alpha channel.
    """
    img = gdal.Open(filename, GA_Update)
    mask = img.GetRasterBand(4)
    data = mask.ReadAsArray(0, 0, mask.XSize, mask.YSize)
    data = np.array([entry * 255 for entry in data])
    mask.WriteArray(data)
    mask.FlushCache()

    # Close the dataset
    img = None
    data = None
    mask = None 

def create_mask_from_fmask(in_l1_dir, out_path):
    log = logging.getLogger(__name__)
    log.info("Creating fmask for {}".format(in_l1_dir))
    with TemporaryDirectory() as td:
        temp_fmask_path = os.path.join(td, "fmask.tif")
        apply_fmask(in_l1_dir, temp_fmask_path)
        fmask_image = gdal.Open(temp_fmask_path)
        fmask_array = fmask_image.GetVirtualMemArray()
        out_image = create_matching_dataset(fmask_image, out_path, datatype=gdal.GDT_Byte)
        out_array = out_image.GetVirtualMemArray(eAccess=gdal.GA_Update)
        log.info("fmask created, converting to binary cloud/shadow mask")
        out_array[:,:] = np.isin(fmask_array, (2, 3, 4), invert=True)
        out_array = None
        out_image = None
        fmask_array = None
        fmask_image = None
        resample_image_in_place(out_path, 10) 

def loadMetadata(self):
        # open the raster and its spatial reference
        self.src = gdal.Open(self.tif_path)

        if self.src is None:
            raise Exception('Could not load GDAL file "%s"' % self.tif_path)
        spatial_reference_raster = osr.SpatialReference(self.src.GetProjection())

        # get the WGS84 spatial reference
        spatial_reference = osr.SpatialReference()
        spatial_reference.ImportFromEPSG(4326)  # WGS84

        # coordinate transformation
        self.coordinate_transform = osr.CoordinateTransformation(spatial_reference, spatial_reference_raster)
        gt = self.geo_transform = self.src.GetGeoTransform()
        dev = (gt[1] * gt[5] - gt[2] * gt[4])
        self.geo_transform_inv = (gt[0], gt[5] / dev, -gt[2] / dev,
                                  gt[3], -gt[4] / dev, gt[1] / dev) 

def get_edge(config, input):

    """ Get sobel edge extraction
    #NOTE: Dependency requirements on the BU Cluster now require this be ran seperately
    #The code below would work (uncommented) if OTB could be loaded
    """

    # The following line creates an instance of the EdgeExtraction application
#    EdgeExtraction = otbApplication.Registry.CreateApplication("EdgeExtraction")

    # The following lines set all the application parameters:
#    EdgeExtraction.SetParameterString("in", input)

#    EdgeExtraction.SetParameterInt("channel", 2) #TODO: set band

    dst_path = config['postprocessing']['deg_class']['edge_dir']
    dst_filename = dst_path + '/' + input.split('/')[-1].split('.')[0] + '_edge.tif'
#    EdgeExtraction.SetParameterString("out", dst_filename)

    # The following line execute the application
#    EdgeExtraction.ExecuteAndWriteOutput()

    im_op = gdal.Open(dst_filename)
    im_ar = im_op.ReadAsArray()
    return im_ar 

def cirrus_correction(stacked_raster_path, out_path):
    stacked_raster = gdal.Open(stacked_raster_path)
    ras_array = stacked_raster.GetVirtualMemArray()
    r_band = ras_array[0]
    g_band = ras_array[1]
    b_band = ras_array[2]
    cirrus_band = ras_array[3]
    out_raster = ras.create_matching_dataset(stacked_raster, out_path, bands = 3)
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)

    for ii, band in enumerate([r_band, g_band, b_band]):
        out_array[ii, ...] = band - ((cirrus_band - 100)*12/(np.log(cirrus_band - 100) +1))

    out_array = None
    b_band = None
    g_band = None
    b_band = None
    cirrus_band = None
    out_array = None
    out_raster = None
    stacked_raster = None 

def count_pixel_classes(map_path, no_data=None):
    """
    Counts pixels in a map. Returns a dictionary of pixels.
    Parameters
    ----------
    map_path: Path to the map to count
    no_data: A value to ignore

    Returns
    -------
    A dictionary of class:count
    """
    map = gdal.Open(map_path)
    map_array = map.GetVirtualMemArray().squeeze()
    unique, counts = np.unique(map_array, return_counts=True)
    out = dict(zip([str(val) for val in unique], counts))
    out.pop(no_data, "_")   # pop the no data value, but don't worry if there's nothing there.
    map_array=None
    map=None
    return out 

def load_file(self, file_path, band_num = 0):
        '''Load the QA file from a give path

        Parameters
            file_path	-	Path of band file.
            band_num	-	Number of band.
        '''
        import gdal

        self.file_path = file_path
        extension = os.path.splitext(file_path)[1].lower()

        # load file according to the file format.
        if extension == '.hdf':
            dataset = gdal.Open(file_path)
            subdataset = dataset.GetSubDatasets()[band_num][0]
            self.band_data = gdal.Open(subdataset).ReadAsArray()
        else:
            bandfile = gdal.Open(file_path)
            self.band_data = bandfile.GetRasterBand(1).ReadAsArray() 

def apply_image_function(in_paths, out_path, function, out_datatype = gdal.GDT_Int32):
    """Applies a pixel-wise function across every image. Assumes each image is exactly contiguous and, for now,
    single-banded. function() should take a list of values and return a single value."""
    rasters = [gdal.Open(in_path) for in_path in in_paths]
    raster_arrays = [raster.GetVirtualMemArray() for raster in rasters]
    in_array = np.stack(raster_arrays, axis=0)

    out_raster = create_matching_dataset(rasters[0], out_path=out_path, datatype=out_datatype)
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    out_array[...] = np.apply_along_axis(function, 0, in_array)

    # Deallocating. Not taking any chances here.
    out_array = None
    out_raster = None
    in_array = None
    for raster_array, raster in zip(raster_arrays, rasters):
        raster_array = None
        raster = None 

def calc_ndvi(raster_path, output_path):
    raster = gdal.Open(raster_path)
    out_raster = create_matching_dataset(raster, output_path, datatype=gdal.GDT_Float32)
    array = raster.GetVirtualMemArray()
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    R = array[2, ...]
    I = array[3, ...]
    out_array[...] = (R-I)/(R+I)

    out_array[...] = np.where(out_array == -2147483648, 0, out_array)

    R = None
    I = None
    array = None
    out_array = None
    raster = None
    out_raster = None 

def raster_to_array(rst_pth):
    """Reads in a raster file and returns a N-dimensional array.

    Parameters
    ----------
    rst_pth
        Path to input raster.
    Returns
    -------
    As N-dimensional array.
    """
    log = logging.getLogger(__name__)
    in_ds = gdal.Open(rst_pth)
    out_array = in_ds.ReadAsArray()

    return out_array 

def get_masked_array(raster, mask_path):
    """
    Returns a numpy.mask masked array for the raster.
    Masked pixels are FALSE in the mask image (multiplicateive map),
    but TRUE in the masked_array (nodata pixels). If the raster is multi-band and
    the mask is single-band, the mask will be applied to every raster.
    Parameters
    ----------
    raster
        A gdal.Image object
    mask_path
        The path to the mask to use

    Returns
    -------
    A numpy.masked array of the raster.

    """
    mask = gdal.Open(mask_path)
    mask_array = mask.GetVirtualMemArray().squeeze()
    raster_array = raster.GetVirtualMemArray()
    # If the shapes do not match, assume single-band mask for multi-band raster
    if len(mask_array.shape) == 2 and len(raster_array.shape) == 3:
        mask_array = project_array(np.asarray(mask_array), raster_array.shape[0], 0)
    return np.ma.array(raster_array, mask=np.logical_not(mask_array)) 

def flatten_probability_image(prob_image, out_path):
    """
    Takes a probability output from classify_image and flattens it into a single layer containing only the maximum
    value from each pixel.

    Parameters
    ----------
    prob_image
        The path to a probability image.
    out_path
        The place to save the flattened image.

    """
    prob_raster = gdal.Open(prob_image)
    out_raster = create_matching_dataset(prob_raster, out_path, bands=1)
    prob_array = prob_raster.GetVirtualMemArray()
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    out_array[:, :] = prob_array.max(axis=0)
    out_array = None
    prob_array = None
    out_raster = None
    prob_raster = None 

def create_mask_from_model(image_path, model_path, model_clear=0, num_chunks=10, buffer_size=0):
    """Returns a multiplicative mask (0 for cloud, shadow or haze, 1 for clear) built from the model at model_path."""
    from pyeo.classification import classify_image  # Deferred import to deal with circular reference
    with TemporaryDirectory() as td:
        log = logging.getLogger(__name__)
        log.info("Building cloud mask for {} with model {}".format(image_path, model_path))
        temp_mask_path = os.path.join(td, "cat_mask.tif")
        classify_image(image_path, model_path, temp_mask_path, num_chunks=num_chunks)
        temp_mask = gdal.Open(temp_mask_path, gdal.GA_Update)
        temp_mask_array = temp_mask.GetVirtualMemArray()
        mask_path = get_mask_path(image_path)
        mask = create_matching_dataset(temp_mask, mask_path, datatype=gdal.GDT_Byte)
        mask_array = mask.GetVirtualMemArray(eAccess=gdal.GF_Write)
        mask_array[:, :] = np.where(temp_mask_array != model_clear, 0, 1)
        temp_mask_array = None
        mask_array = None
        temp_mask = None
        mask = None
        if buffer_size:
            buffer_mask_in_place(mask_path, buffer_size)
        log.info("Cloud mask for {} saved in {}".format(image_path, mask_path))
        return mask_path 

def strip_bands(in_raster_path, out_raster_path, bands_to_strip):
    in_raster = gdal.Open(in_raster_path)
    out_raster_band_count = in_raster.RasterCount-len(bands_to_strip)
    out_raster = create_matching_dataset(in_raster, out_raster_path, bands=out_raster_band_count)
    out_raster_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    in_raster_array = in_raster.GetVirtualMemArray()

    bands_to_copy = [band for band in range(in_raster_array.shape[0]) if band not in bands_to_strip]

    out_raster_array[...] = in_raster_array[bands_to_copy, :,:]

    out_raster_array = None
    in_raster_array = None
    out_raster = None
    in_raster = None

    return out_raster_path 

def create_mask_from_class_map(class_map_path, out_path, classes_of_interest, buffer_size=0, out_resolution=None):
    """Creates a mask from a classification mask: 1 for each pixel containing one of classes_of_interest, otherwise 0"""
    # TODO: pull this out of the above function
    class_image = gdal.Open(class_map_path)
    class_array = class_image.GetVirtualMemArray()
    mask_array = np.isin(class_array, classes_of_interest)
    out_mask = create_matching_dataset(class_image, out_path, datatype=gdal.GDT_Byte)
    out_array = out_mask.GetVirtualMemArray(eAccess=gdal.GA_Update)
    np.copyto(out_array, mask_array)
    class_array = None
    class_image = None
    out_array = None
    out_mask = None
    if out_resolution:
        resample_image_in_place(out_path, out_resolution)
    if buffer_size:
        buffer_mask_in_place(out_path, buffer_size)
    return out_path 

def test_pixel_bounds_from_polygon(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    raster = gdal.Open(os.path.join(test_dir.path, "ne_test"))
    aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    layer = aoi.GetLayer(0)
    aoi_feature = layer.GetFeature(0)
    polygon = aoi_feature.GetGeometryRef()
    result = pyeo.coordinate_manipulation.pixel_bounds_from_polygon(raster, polygon)
    assert result == (1, 11, 1, 11) 

def test_point_to_pixel_coordinates(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    point = r"POINT (101 112)"
    raster = gdal.Open(os.path.join(test_dir.path, "ne_test"))
    out = pyeo.coordinate_manipulation.point_to_pixel_coordinates(raster, point)
    # Since ne_test is an 11 by 12 image, tl corner 0,0 w 10m pixels,
    # we'd expect the br coords to be 10, 11
    assert out == (10, 11) 

def buffer_mask_in_place(mask_path, buffer_size):
    """Expands a mask in-place, overwriting the previous mask"""
    log = logging.getLogger(__name__)
    log.info("Buffering {} with buffer size {}".format(mask_path, buffer_size))
    mask = gdal.Open(mask_path, gdal.GA_Update)
    mask_array = mask.GetVirtualMemArray(eAccess=gdal.GA_Update)
    cache = morph.binary_erosion(mask_array.squeeze(), selem=morph.disk(buffer_size))
    np.copyto(mask_array, cache)
    mask_array = None
    mask = None 

def test_reprojection():
    os.chdir(os.path.dirname(os.path.abspath(__file__)))
    try:
        os.remove(r"test_outputs/reprojection_test.tif")
    except FileNotFoundError:
        pass
    new_projection = r"""PROJCS["WGS 84 / UTM zone 36S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",33],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32736"]]"""
    image = r"test_data/composite_T36MZE_20190509T073621_20190519T073621_clipped.tif"
    out_file = r"test_outputs/reprojection_test.tif"
    pyeo.raster_manipulation.reproject_image(image, out_file, new_projection, do_post_resample=False)
    result = gdal.Open(out_file)
    assert result
   # assert result.GetProjection() == new_projection
    result_array = result.GetVirtualMemArray()
    assert result_array.max() > 10 

def test_clip_raster(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_raster = os.path.join(test_dir.path, "se_test")
    test_aoi = os.path.join(test_dir.path, "aoi")
    result_path = os.path.join(test_dir.path, "test_out.tif")
    pyeo.raster_manipulation.clip_raster(test_raster, test_aoi, result_path)
    result = gdal.Open(result_path)
    assert result.ReadAsArray().shape == (5, 10, 10)
    assert result.GetGeoTransform() == (10, 10, 0, 10, 0, -10) 

def test_check_overlap(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_raster = gdal.Open(os.path.join(test_dir.path, "ne_test"))
    test_aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    assert pyeo.coordinate_manipulation.check_overlap(test_raster, test_aoi) is True 

def test_get_aoi_intersection(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_raster = gdal.Open(os.path.join(test_dir.path, "se_test"))
    test_aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    result = pyeo.coordinate_manipulation.get_aoi_intersection(test_raster, test_aoi)
    assert result.GetEnvelope() == (10, 20, 10, 20) 

def test_reshape_raster_for_ml(managed_ml_geotiff_dir):
    test_dir = managed_ml_geotiff_dir
    test_image_path = os.path.join(test_dir.path, "training_image")
    test_image = gdal.Open(test_image_path)
    array = test_image.GetVirtualMemArray()
    array = pyeo.classification.reshape_raster_for_ml(array)
    assert array.shape == (120, 8)
    assert np.all(array[2, :] == [2]*8)
    assert np.all(array[119, :] == [119]*8) 

def test_aoi_to_mask(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    result_path = os.path.join(test_dir.path, "test_out.tif")
    pyeo.aoi_to_mask(test_aoi, result_path)
    result = gdal.Open(result_path)
    assert result.ReadAsArray()[2, 2] == 1