Python osgeo.gdal.Open() Examples

def GetGeoInfo(FileName):

    if exists(FileName) is False:
            raise Exception('[Errno 2] No such file or directory: \'' + FileName + '\'')    
    
    
    SourceDS = gdal.Open(FileName, gdal.GA_ReadOnly)
    if SourceDS == None:
        raise Exception("Unable to read the data file")
    
    NDV = SourceDS.GetRasterBand(1).GetNoDataValue()
    xsize = SourceDS.RasterXSize
    ysize = SourceDS.RasterYSize
    GeoT = SourceDS.GetGeoTransform()
    Projection = osr.SpatialReference()
    Projection.ImportFromWkt(SourceDS.GetProjectionRef())
    DataType = SourceDS.GetRasterBand(1).DataType
    DataType = gdal.GetDataTypeName(DataType)
    
    return NDV, xsize, ysize, GeoT, Projection, DataType
#==============================================================================

#==============================================================================
# Function to read the original file's projection: 

def as_array(path, band_axis=True):
    '''
    A convenience function for opening a raster as an array and accessing its
    spatial information; takes a single string argument. Arguments:
        path        The path of the raster file to open as an array
        band_axis   True to include a band axis, even for single-band rasters
    '''
    ds = gdal.Open(path)
    arr = ds.ReadAsArray()
    gt = ds.GetGeoTransform()
    wkt = ds.GetProjection()
    ds = None

    # Make sure that single-band rasters have a band axis
    if band_axis and len(arr.shape) < 3:
        shp = arr.shape
        arr = arr.reshape((1, shp[0], shp[1]))

    return (arr, gt, wkt) 

def fix_pyramid(scene_root, files, verbose=False):
    ret = False
    for filename in files:
        if not filename.endswith('.TIF'):
            continue

        ds = gdal.Open(filename)
        ov_count = ds.GetRasterBand(1).GetOverviewCount()
        ds = None

        if ov_count != 0:
            continue
        
        ret = True
        
        splitter.build_pyramid(filename, verbose=verbose)

    return ret 

def fix_tiling(scene_root, files, verbose=False):
    ret = False
    for filename in files:
        if not filename.endswith('.TIF'):
            continue

        ds = gdal.Open(filename)
        bx, by = ds.GetRasterBand(1).GetBlockSize()
        ds = None
        
        if by != 1:
            continue

        ret = True
        
        splitter.internally_compress(filename, verbose=verbose)

    return ret 

def loadProductOptical(filename):
    import numpy as np
    '''
    Load the product using Product Manager.
    '''
    ds = gdal.Open(filename)
#    pdb.set_trace()
    band = ds.GetRasterBand(1)
    
    img = band.ReadAsArray()
    img = img.astype(np.float32)
    
    band=None
    ds=None
    
    return img 

def test_array_to_raster_interface(self):
        '''
        The array_to_raster() and array_to_raster_clone functions should
        perform as expected.
        '''
        # First round
        ds = gdal.Open(os.path.join(self.test_dir, 'multi3_raster.tiff'))
        gt = ds.GetGeoTransform()
        wkt = ds.GetProjection()
        arr = ds.ReadAsArray()
        ds = None
        rast = array_to_raster(arr, gt, wkt)
        self.assertEqual(gt, rast.GetGeoTransform())
        self.assertEqual(wkt, rast.GetProjection())

        # Second round
        rast2 = array_to_raster_clone(arr, os.path.join(self.test_dir,
            'multi7_raster.tiff'))
        self.assertEqual(gt, rast2.GetGeoTransform())
        self.assertEqual(wkt, rast2.GetProjection()) 

def __init__(self, sentinel1_product_zip_path):
        self.product_zip_path = sentinel1_product_zip_path
        sentinel1_product_dir = os.path.dirname(sentinel1_product_zip_path)
        sentinel1_product_zipname = os.path.basename(sentinel1_product_zip_path)
        self.product_dir = sentinel1_product_dir
        self.sentinel1_product_zipname = sentinel1_product_zipname
        self.granule_identifier, _ = os.path.splitext(sentinel1_product_zipname)

        manifest_path = extract_manifest_from_zip(sentinel1_product_zip_path)
        self.manifest_tree = ET.parse(manifest_path)
        try:
            os.remove(manifest_path)
            os.rmdir(os.path.dirname(manifest_path))
        except:
            log.warn("Cannot cleanup manifest directory")

        #sentinel1_safe_pkg_path = "/vsizip/{}/{}.SAFE/manifest.safe".format(self.product_zip_path, self.granule_identifier)

        #self.safe_package_path = sentinel1_safe_pkg_path
        #self.datastore = gdal.Open(sentinel1_safe_pkg_path) 

def loadProductOptical(filename):
    import numpy as np
    '''
    Load the product using Product Manager.
    '''
    ds = gdal.Open(filename)
#    pdb.set_trace()
    band = ds.GetRasterBand(1)
    
    img = band.ReadAsArray()
    img = img.astype(np.float32)
    
    band=None
    ds=None
    
    return img 

def get_metadata(self):
        manifest_zip_path = get_manifest_zip_path(self.product_zip_path)
        datastore = gdal.Open(manifest_zip_path)
        metadata_dict = datastore.GetMetadata()
        metadata_dict['NAME'] = self.granule_identifier
        startTime = metadata_dict['ACQUISITION_START_TIME']
        endTime   = metadata_dict['ACQUISITION_STOP_TIME']
        # round to milliseconds
        m = re.search("(.*)(\d{3})(\d{3})$", startTime)
        if m:
            startTime = m.groups()[0] + m.groups()[1]
        m = re.search("(.*)(\d{3})(\d{3})$", endTime)
        if m:
            endTime = m.groups()[0] + m.groups()[1]
        metadata_dict['ACQUISITION_START_TIME'] = startTime + 'Z'
        metadata_dict['ACQUISITION_STOP_TIME'] = endTime + 'Z'
        return metadata_dict 

def read_image():
    """ Read image ``f`` and return ``np.array`` of image values

    Image will be (nband, nrow, ncol)
    """
    def _read_image(f):
        ds = gdal.Open(f, gdal.GA_ReadOnly)
        dtype = gdal_array.GDALTypeCodeToNumericTypeCode(
            ds.GetRasterBand(1).DataType)
        nrow, ncol, nband = ds.RasterYSize, ds.RasterYSize, ds.RasterCount
        img = np.empty((nband, nrow, ncol), dtype=dtype)
        for i_b in range(nband):
            b = ds.GetRasterBand(i_b + 1)
            img[i_b, ...] = b.ReadAsArray()
        return img
    return _read_image 

def test_spectral_profile(self):
        '''
        Should correctly retrieve a spectral profile from a raster dataset.
        '''
        coords = ((-84.5983, 42.7256), (-85.0807, 41.1138))
        pixels = [(18, 0), (2, 59)]
        file_path = os.path.join(self.test_dir, 'multi3_raster.tiff')
        ds = gdal.Open(file_path)
        kwargs = {
            'gt': ds.GetGeoTransform(),
            'wkt': ds.GetProjection(),
            'dd': True
        }

        # The true spectral profile
        spectra = np.array([[237, 418, 325], [507, 616, 445]], dtype=np.int16)
        sp1 = spectra_at_xy(ds, coords, **kwargs)
        sp2 = spectra_at_xy(ds.ReadAsArray(), coords, **kwargs)
        sp3 = spectra_at_idx(ds.ReadAsArray().transpose(), pixels)
        self.assertEqual(spectra.tolist(), sp1.tolist())
        self.assertEqual(spectra.tolist(), sp2.tolist())
        self.assertEqual(spectra.tolist(), sp3.tolist()) 

def get_image_attribute(image_filename):
    """ Use GDAL to open image and return some attributes

    Args:
        image_filename (str): image filename

    Returns:
        tuple: nrow (int), ncol (int), nband (int), NumPy datatype (type)
    """
    try:
        image_ds = gdal.Open(image_filename, gdal.GA_ReadOnly)
    except Exception as e:
        logger.error('Could not open example image dataset ({f}): {e}'
                     .format(f=image_filename, e=str(e)))
        raise

    nrow = image_ds.RasterYSize
    ncol = image_ds.RasterXSize
    nband = image_ds.RasterCount
    dtype = gdal_array.GDALTypeCodeToNumericTypeCode(
        image_ds.GetRasterBand(1).DataType)

    return (nrow, ncol, nband, dtype) 

def test_hall_rectification(self):
        '''Should rectify an image in the expected way.'''
        control_sets = {
            'High/Bright': [(331501.45,4694346.66), (319495.39,4706820.66), (298527.006,4691417.99)],
            'Low/Dark': [(322577.40,4658508.99), (361612.79,4694665.62), (378823.69,4692132.56)]
        }
        ref_raster = gdal.Open(os.path.join(self.test_dir, 'multi7_raster.tiff'))
        sub_raster = gdal.Open(os.path.join(self.test_dir, 'multi7_raster2.tiff'))

        # NOTE: Using same control sets for reference, subject
        hall_rectification(ref_raster, sub_raster, self.test_dir, control_sets, control_sets)

        arr, gt, wkt = as_array(os.path.join(self.test_dir, 'rect_multi7_raster2.tiff'))
        self.assertTrue(np.array_equal(arr.shape, (6, 74, 81)))
        self.assertTrue(np.array_equiv(arr[:,50,50].round(5), np.array([
            17, 1331, 1442, 3422, 2916, 2708
        ]).round(5))) 

def chunk(raster_filename, chunk_size=256, chunk_dir='/tmp/'):
    """
    Given a raster, a chunk size, and a directory to write into...
    Break the raster up into chunks of the appropriate size.
    """
    CROP_CMD = 'gdal_translate -co ALPHA=YES -co PHOTOMETRIC=RGB -srcwin %s %s %s %s %s %s'
    # % (xoff, yoff, xsize, ysize, src, dst)

    base = os.path.basename(os.path.splitext(raster_filename)[0])

    ds = gdal.Open(raster_filename)
    numPixelsWide, numPixelsHigh = ds.RasterXSize, ds.RasterYSize
    for x in range(0, numPixelsWide-chunk_size-1, chunk_size):
        for y in range(0, numPixelsHigh-chunk_size-1, chunk_size):
            chunk_filename = os.path.join(
                chunk_dir, '%s-%s-%s.tif' % (base, x, y)
            )
            os.system(CROP_CMD % (
                x, y, chunk_size, chunk_size, raster_filename, chunk_filename
            ))
            yield chunk_filename 

def getNoDataValue(rasterfn):
    """This gets the nodata value from the raster

    Args:
        rasterfn (str): The filename (with path and extension) of the raster

    Returns:
        float: nodatavalue; the nodata value

    Author: SMM
    """
    raster = gdal.Open(rasterfn)
    band = raster.GetRasterBand(1)
    return band.GetNoDataValue()
#==============================================================================

#============================================================================== 

def setNoDataValue(rasterfn):
    """This sets the nodata value from the raster

    Args:
        rasterfn (str): The filename (with path and extension) of the raster

    Returns:
        None

    Author: SMM
    """

    raster = gdal.Open(rasterfn)
    band = raster.GetRasterBand(1)
    return band.SetNoDataValue()
#==============================================================================

#============================================================================== 

def get_extent(srcFileImage):
    gdata = gdal.Open(srcFileImage)
    geo = gdata.GetGeoTransform()
    # data = gdata.ReadAsArray()

    xres = geo[1]
    yres = geo[5]
    # xmin = geo[0]
    # xmax = geo[0] + (xres * gdata.RasterXSize)
    # ymin = geo[3] + (yres * gdata.RasterYSize)
    # ymax = geo[3]
    xmin = geo[0] + xres * 0.5
    xmax = geo[0] + (xres * gdata.RasterXSize) - xres * 0.5
    ymin = geo[3] + (yres * gdata.RasterYSize) + yres * 0.5
    ymax = geo[3] - yres * 0.5

    return xmin, ymin, xmax, ymax


############################################################################### 

def get_gsd(im_test_file):
    '''return gsd in meters'''
    srcImage = gdal.Open(im_test_file)
    geoTrans = srcImage.GetGeoTransform()
    ulX = geoTrans[0]
    ulY = geoTrans[3]
    # xDist = geoTrans[1]
    yDist = geoTrans[5]
    # rtnX = geoTrans[2]
    # rtnY = geoTrans[4]

    # get haversine distance
    # dx = _haversine(ulX, ulY, ulX+xDist, ulY) #haversine(lon1, lat1, lon2, lat2)
    dy = _haversine(ulX, ulY, ulX, ulY+yDist)   #haversine(lon1, lat1, lon2, lat2)

    return dy  # dx


############################################################################### 

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 

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 

def create_mask_from_vector(vector_data_path, cols, rows, geo_transform, 
			    projection, target_value=1,
                            output_fname='', dataset_format='MEM'):

    """
    Rasterize the given vector (wrapper for gdal.RasterizeLayer). 
    Return a gdal.Dataset.
    :param vector_data_path: Path to a shapefile
    :param cols: Number of columns of the result
    :param rows: Number of rows of the result
    :param geo_transform: Returned value of gdal.Dataset.GetGeoTransform 
	(coefficients for transforming between pixel/line (P,L) raster space,
	 and projection coordinates (Xp,Yp) space.
    :param projection: Projection definition string (Returned by 
	gdal.Dataset.GetProjectionRef)
    :param target_value: Pixel value for the pixels. Must be a valid 
	gdal.GDT_UInt16 value.
    :param output_fname: If the dataset_format is GeoTIFF, this is the output 
	file name
    :param dataset_format: The gdal.Dataset driver name. [default: MEM]
    """

    driver = ogr.GetDriverByName('ESRI Shapefile')
    data_source = driver.Open(vector_data_path, 0)
    if data_source is None:
        report_and_exit("File read failed: %s", vector_data_path)
    layer = data_source.GetLayer(0)
    driver = gdal.GetDriverByName(dataset_format)
    target_ds = driver.Create(output_fname, cols, rows, 1, gdal.GDT_UInt16)
    target_ds.SetGeoTransform(geo_transform)
    target_ds.SetProjection(projection)
    gdal.RasterizeLayer(target_ds, [1], layer, burn_values=[target_value])
    return target_ds 

def get_img_array_from_zip(zip_file, img_name):
  with tempfile.NamedTemporaryFile() as temp_file:
    # Copy to temp file.
    with zip_file.open(img_name) as img_file:
        shutil.copyfileobj(img_file, temp_file)

    # Extract as numpy array.
    geo = gdal.Open(temp_file.name)
    return geo.ReadAsArray() if geo is not None else None 

def test_get_pixel_latlon():
    # Expected out for  top-left corner of test image, (0,0)
    # Test image is in EPSG 32748, QGIS says that TL corner coords are 600001.8, 9399997.9
    # epsg.io says that this is 105.9026743, -5.4275703 in latlon
    os.chdir(pathlib.Path(__file__).parent)
    test_image_path = "test_data/S2A_MSIL2A_20170922T025541_N0205_R032_T48MXU_20170922T031450.SAFE/GRANULE/L2A_T48MXU_A011755_20170922T031450/IMG_DATA/R20m/L2A_T48MXU_20170922T025541_AOT_20m.jp2"
    target_lon = 105.9026743
    target_lat = -5.4275703
    test_image = gdal.Open(test_image_path)
    out_lat, out_lon = terrain_correction.get_pixel_latlon(test_image, 0, 0)
    np.testing.assert_allclose(out_lat, target_lat, 0.001)
    np.testing.assert_allclose(out_lon, target_lon, 0.001) 

def test_terrain_correction_landsat(monkeypatch):
    dem_path = "test_data/dem_test_indonesia.tif"
    in_path = "test_data/landsat_stack.tif"
    raster_timezone = pytz.timezone("UTC")
    raster_datetime = dt.datetime(2015, 7, 5, 3, 5, 42, tzinfo=raster_timezone)
    out_path = "test_outputs/correction_landsat_indonesia.tif"
    terrain_correction.do_terrain_correction(in_path, dem_path, out_path, raster_datetime, is_landsat=False)
    assert gdal.Open(out_path) 

def test_get_dem_slope_and_angle():
    test_dem_path = r"test_data/N001E109/N001E109_AVE_DSM.tif"
    slope_path = r"test_outputs/slope.dem"
    angle_path = r"test_outputs/angle.dem"
    terrain_correction.get_dem_slope_and_angle(test_dem_path, slope_path, angle_path)
    assert gdal.Open(slope_path)
    assert gdal.Open(angle_path) 

def test_calculate_latlon_array():
    raster_path = "test_data/dem_test_indonesia.tif"
    raster = gdal.Open(raster_path)
    array = raster.GetVirtualMemArray()
    transformer, gt = terrain_correction._generate_latlon_transformer(raster)
    lon, lat = terrain_correction._generate_latlon_arrays(array, transformer, gt)
    test_lat = joblib.load("test_data/lat_array_indo")
    test_lon = joblib.load("test_data/lon_array_indo")
    assert np.all(lat == test_lat)
    assert np.all(lon == test_lon) 

def test_masking(self):
        '''
        Masking should go on without a hitch and the result should be just
        as expected.
        '''
        file_path = os.path.join(self.test_dir, 'multi7_raster.tiff')
        ds = gdal.Open(file_path)
        raw_mask = gdal.Open(os.path.join(self.test_dir, 'multi7_mask.tiff'))
        mask = cfmask(raw_mask, nodata=-9999)
        masked = binary_mask(ds.ReadAsArray(), mask)
        self.assertEqual(ds.ReadAsArray().diagonal()[0,0], 0)
        self.assertEqual(masked.diagonal()[0,0], -9999) 

def _latlon2pixel(lat, lon, input_raster='', targetsr='', geom_transform=''):
    '''
    Convert latitude, longitude coords to pixexl coords.
    From spacenet geotools
    '''

    sourcesr = osr.SpatialReference()
    sourcesr.ImportFromEPSG(4326)

    geom = ogr.Geometry(ogr.wkbPoint)
    # geom.AddPoint(lon, lat)
    geom.AddPoint(lat, lon)

    if targetsr == '':
        src_raster = gdal.Open(input_raster)
        targetsr = osr.SpatialReference()
        targetsr.ImportFromWkt(src_raster.GetProjectionRef())
    coord_trans = osr.CoordinateTransformation(sourcesr, targetsr)
    if geom_transform == '':
        src_raster = gdal.Open(input_raster)
        transform = src_raster.GetGeoTransform()
    else:
        transform = geom_transform

    x_origin = transform[0]
    # print(x_origin)
    y_origin = transform[3]
    # print(y_origin)
    pixel_width = transform[1]
    # print(pixel_width)
    pixel_height = transform[5]
    # print(pixel_height)
    geom.Transform(coord_trans)
    # print(geom.GetPoint())
    x_pix = (geom.GetPoint()[0] - x_origin) / pixel_width
    y_pix = (geom.GetPoint()[1] - y_origin) / pixel_height

    return (x_pix, y_pix)


############################################################################### 

def _getLOS(self, filename, component):
        path = op.dirname(filename)
        fn = glob.glob(op.join(path, component))
        if len(fn) != 1:
            raise ImportError('Cannot find LOS vector file %s!' % component)

        dataset = gdal.Open(fn[0], gdal.GA_ReadOnly)
        return self._readBandData(dataset) 

def test_terrain_correction_s2():
    dem_path = "test_data/dem_test_clipped.tif"
    in_path = "test_data/S2A_MSIL1C_20200310T025541_N0209_R032_T48MXU_20200310T062815.tif"
    out_path = "test_outputs/correction_s2_indonesia.tif"
    if os.path.exists(out_path):
        os.remove(out_path)
    raster_timezone = pytz.timezone("UTC")
    raster_datetime = dt.datetime(2020, 3, 10, 2, 55, 41, tzinfo=raster_timezone)
    terrain_correction.do_terrain_correction(in_path, dem_path, out_path, raster_datetime)
    out = gdal.Open(out_path)
    assert out.GetVirtualMemArray().max() > 10