Python gdal.Open方法代码示例

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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() 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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() 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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)) 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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) 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Open [as 别名]
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