Python osgeo.gdal.Open() Examples

The following are 30 code examples of osgeo.gdal.Open(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module osgeo.gdal , or try the search function .
Example #1
Source Project: LSDMappingTools   Author: LSDtopotools   File: LSDMappingTools.py    License: MIT License 7 votes vote down vote up
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: 
Example #2
Source Project: yatsm   Author: ceholden   File: conftest.py    License: MIT License 6 votes vote down vote up
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 
Example #3
Source Project: yatsm   Author: ceholden   File: readers.py    License: MIT License 6 votes vote down vote up
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) 
Example #4
Source Project: cloudless   Author: BradNeuberg   File: populate_db.py    License: Apache License 2.0 6 votes vote down vote up
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 
Example #5
Source Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_GDALIO.py    License: MIT License 6 votes vote down vote up
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()
#==============================================================================

#============================================================================== 
Example #6
Source Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_GDALIO.py    License: MIT License 6 votes vote down vote up
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()
#==============================================================================

#============================================================================== 
Example #7
Source Project: lidar   Author: giswqs   File: slicing.py    License: 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 #8
Source Project: unmixing   Author: arthur-e   File: utils.py    License: MIT License 6 votes vote down vote up
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) 
Example #9
Source Project: unmixing   Author: arthur-e   File: tests.py    License: MIT License 6 votes vote down vote up
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))) 
Example #10
Source Project: unmixing   Author: arthur-e   File: tests.py    License: MIT License 6 votes vote down vote up
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()) 
Example #11
Source Project: unmixing   Author: arthur-e   File: tests.py    License: MIT License 6 votes vote down vote up
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()) 
Example #12
Source Project: apls   Author: CosmiQ   File: apls_utils.py    License: Apache License 2.0 6 votes vote down vote up
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


############################################################################### 
Example #13
Source Project: apls   Author: CosmiQ   File: apls_utils.py    License: Apache License 2.0 6 votes vote down vote up
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


############################################################################### 
Example #14
Source Project: autoRIFT   Author: leiyangleon   File: testautoRIFT_ISCE.py    License: Apache License 2.0 6 votes vote down vote up
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 
Example #15
Source Project: autoRIFT   Author: leiyangleon   File: testautoRIFT.py    License: Apache License 2.0 6 votes vote down vote up
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 
Example #16
Source Project: landsat_ingestor   Author: landsat-pds   File: reprocess_scene.py    License: Apache License 2.0 6 votes vote down vote up
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 
Example #17
Source Project: landsat_ingestor   Author: landsat-pds   File: reprocess_scene.py    License: Apache License 2.0 6 votes vote down vote up
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 
Example #18
Source Project: evo-odas   Author: geosolutions-it   File: S1Reader.py    License: MIT License 6 votes vote down vote up
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) 
Example #19
Source Project: evo-odas   Author: geosolutions-it   File: S1Reader.py    License: MIT License 6 votes vote down vote up
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 
Example #20
Source Project: fetchLandsatSentinelFromGoogleCloud   Author: vascobnunes   File: fels.py    License: MIT License 5 votes vote down vote up
def check_full_tile(image):
    gdalData = gdal.Open(image)
    if gdalData is None:
        sys.exit("ERROR: can't open raster")

    # get width and heights of the raster
    xsize = gdalData.RasterXSize
    ysize = gdalData.RasterYSize

    # process the raster
    band_i = gdalData.GetRasterBand(1)
    raster = band_i.ReadAsArray()

    # create dictionary for unique values count
    count = {}

    # count unique values for the given band
    for col in range(xsize):
        for row in range(ysize):
            cell_value = raster[row, col]

            # check if cell_value is NaN
            if cell_value == 0:
                # add cell_value to dictionary
                if cell_value in count:
                    count[cell_value] += 1
                else:
                    count[cell_value] = 1
                break
    for key in sorted(count.keys()):
        if count[key] is not None:
            return "Partial" 
Example #21
Source Project: radiometric_normalization   Author: planetlabs   File: gimage.py    License: Apache License 2.0 5 votes vote down vote up
def load(filename, nodata=None, last_band_alpha=False):
    logging.info('GImage: Loading {} as GImage'.format(filename))
    gdal_ds = gdal.Open(filename)
    if gdal_ds is None:
        raise Exception('Unable to open file "{}" with gdal.Open()'.format(
            filename))

    alpha, band_count = read_alpha_and_band_count(gdal_ds, last_band_alpha)
    bands = _read_all_bands(gdal_ds, band_count)
    metadata = read_metadata(gdal_ds)

    if nodata:
        alpha = alpha * _nodata_to_mask(bands, nodata)
    return GImage(bands, alpha, metadata) 
Example #22
Source Project: radiometric_normalization   Author: planetlabs   File: normalize_wrapper.py    License: Apache License 2.0 5 votes vote down vote up
def _open_image_and_get_info(path, last_band_alpha):
    gdal_ds = gdal.Open(path)
    alpha_band, band_count = gimage.read_alpha_and_band_count(
        gdal_ds, last_band_alpha=last_band_alpha)
    return gdal_ds, alpha_band, band_count 
Example #23
Source Project: radiometric_normalization   Author: planetlabs   File: display_wrapper.py    License: Apache License 2.0 5 votes vote down vote up
def _open_image_and_get_info(path, last_band_alpha):
    gdal_ds = gdal.Open(path)
    alpha_band, band_count = gimage.read_alpha_and_band_count(
        gdal_ds, last_band_alpha=last_band_alpha)
    return gdal_ds, alpha_band, band_count 
Example #24
Source Project: radiometric_normalization   Author: planetlabs   File: transformation_wrapper.py    License: Apache License 2.0 5 votes vote down vote up
def _open_image_and_get_info(path, last_band_alpha):
    gdal_ds = gdal.Open(path)
    alpha_band, band_count = gimage.read_alpha_and_band_count(
        gdal_ds, last_band_alpha=last_band_alpha)
    return gdal_ds, alpha_band, band_count 
Example #25
Source Project: radiometric_normalization   Author: planetlabs   File: pif_wrapper.py    License: Apache License 2.0 5 votes vote down vote up
def _open_image_and_get_info(path, last_band_alpha):
    gdal_ds = gdal.Open(path)
    alpha_band, band_count = gimage.read_alpha_and_band_count(
        gdal_ds, last_band_alpha=last_band_alpha)
    return gdal_ds, alpha_band, band_count 
Example #26
Source Project: radiometric_normalization   Author: planetlabs   File: gimage_tests.py    License: Apache License 2.0 5 votes vote down vote up
def test__read_all_bands(self):
        gdal_ds = gdal.Open(self.test_photometric_alpha_image)
        bands = gimage._read_all_bands(gdal_ds, 3)
        numpy.testing.assert_array_equal(bands[0], self.band)
        numpy.testing.assert_array_equal(bands[1], self.band)
        numpy.testing.assert_array_equal(bands[2], self.band) 
Example #27
Source Project: radiometric_normalization   Author: planetlabs   File: gimage_tests.py    License: Apache License 2.0 5 votes vote down vote up
def test_read_single_band(self):
        gdal_ds = gdal.Open(self.test_photometric_alpha_image)
        band = gimage.read_single_band(gdal_ds, 3)
        numpy.testing.assert_array_equal(band, self.band) 
Example #28
Source Project: radiometric_normalization   Author: planetlabs   File: gimage_tests.py    License: Apache License 2.0 5 votes vote down vote up
def test_read_alpha_and_band_count(self):
        gdal_ds = gdal.Open(self.test_photometric_alpha_image)
        alpha, band_count = gimage.read_alpha_and_band_count(gdal_ds)

        self.assertEqual(band_count, 3)
        numpy.testing.assert_array_equal(alpha, self.mask) 
Example #29
Source Project: radiometric_normalization   Author: planetlabs   File: gimage_tests.py    License: Apache License 2.0 5 votes vote down vote up
def test__save_to_ds(self):
        output_file = 'test_save_to_ds.tif'

        test_band = numpy.array([[0, 1], [2, 3]], dtype=numpy.uint16)
        test_gimage = gimage.GImage([test_band], self.mask, self.metadata)
        output_ds = gdal.GetDriverByName('GTiff').Create(
            output_file, 2, 2, 2, gdal.GDT_UInt16,
            options=['ALPHA=YES'])
        gimage._save_to_ds(test_gimage, output_ds, nodata=3)

        # Required for gdal to write to file
        output_ds = None

        test_ds = gdal.Open(output_file)

        saved_number_of_bands = test_ds.RasterCount
        self.assertEquals(saved_number_of_bands, 2)

        saved_band = test_ds.GetRasterBand(1).ReadAsArray()
        numpy.testing.assert_array_equal(saved_band, self.band)

        saved_nodata = test_ds.GetRasterBand(1).GetNoDataValue()
        self.assertEqual(saved_nodata, 3)

        saved_alpha = test_ds.GetRasterBand(2).ReadAsArray()
        numpy.testing.assert_array_equal(saved_alpha, self.mask * 255)

        os.unlink(output_file) 
Example #30
Source Project: yatsm   Author: ceholden   File: readers.py    License: MIT License 5 votes vote down vote up
def read_pixel_timeseries(files, px, py):
    """ Returns NumPy array containing timeseries values for one pixel

    Args:
        files (list): List of filenames to read from
        px (int): Pixel X location
        py (int): Pixel Y location

    Returns:
        np.ndarray: Array (nband x n_images) containing all timeseries data
            from one pixel
    """
    nrow, ncol, nband, dtype = get_image_attribute(files[0])

    if px < 0 or px >= ncol or py < 0 or py >= nrow:
        raise IndexError('Row/column {r}/{c} is outside of image '
                         '(nrow/ncol: {nrow}/{ncol})'.
                         format(r=py, c=px, nrow=nrow, ncol=ncol))

    Y = np.zeros((nband, len(files)), dtype=dtype)

    for i, f in enumerate(files):
        ds = gdal.Open(f, gdal.GA_ReadOnly)
        for b in range(nband):
            Y[b, i] = ds.GetRasterBand(b + 1).ReadAsArray(px, py, 1, 1)

    return Y