Python gdal.Open方法代碼示例

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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: 

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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 

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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) 

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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 

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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()
#==============================================================================

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

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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()
#==============================================================================

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

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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 

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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) 

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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))) 

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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()) 

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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()) 

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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


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

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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


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

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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 

# 需要導入模塊: from osgeo import gdal [as 別名]
# 或者: from osgeo.gdal import Open [as 別名]
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