Python osgeo.gdal 模块，GA_ReadOnly() 实例源码

def read_naip(file_path, bands_to_use):
    """
    Read in a NAIP, based on www.machinalis.com/blog/python-for-geospatial-data-processing.

    Bands_to_use is an array like [0,0,0,1], designating whether to use each band (R, G, B, IR).
    """
    raster_dataset = gdal.Open(file_path, gdal.GA_ReadOnly)

    bands_data = []
    index = 0
    for b in range(1, raster_dataset.RasterCount + 1):
        band = raster_dataset.GetRasterBand(b)
        if bands_to_use[index] == 1:
            bands_data.append(band.ReadAsArray())
        index += 1
    bands_data = numpy.dstack(bands_data)

    return raster_dataset, bands_data

def tag_with_locations(test_images, predictions, tile_size, state_abbrev):
    """Combine image data with label data, so info can be rendered in a map and list UI.

    Add location data for convenience too.
    """
    combined_data = []
    for idx, img_loc_tuple in enumerate(test_images):
        raster_filename = img_loc_tuple[2]
        raster_dataset = gdal.Open(os.path.join(NAIP_DATA_DIR, raster_filename), gdal.GA_ReadOnly)
        raster_tile_x = img_loc_tuple[1][0]
        raster_tile_y = img_loc_tuple[1][1]
        ne_lat, ne_lon = pixel_to_lat_lon_web_mercator(raster_dataset, raster_tile_x +
                                                       tile_size, raster_tile_y)
        sw_lat, sw_lon = pixel_to_lat_lon_web_mercator(raster_dataset, raster_tile_x,
                                                       raster_tile_y + tile_size)
        certainty = predictions[idx][0]
        formatted_info = {'certainty': certainty, 'ne_lat': ne_lat, 'ne_lon': ne_lon,
                          'sw_lat': sw_lat, 'sw_lon': sw_lon, 'raster_tile_x': raster_tile_x,
                          'raster_tile_y': raster_tile_y, 'raster_filename': raster_filename,
                          'state_abbrev': state_abbrev, 'country_abbrev': 'USA'
                          }
        combined_data.append(formatted_info)
    return combined_data

def readImageMetaData(inputFile, name):
    # Open the dataset in read only mode
    dataset = gdal.Open(inputFile, gdal.GA_ReadOnly)
    # Check that the image has been opened
    if not dataset is None:
        # Get the meta value specified
        metaData = dataset.GetMetaDataItem(name)
        # Check that it is present
        if metaData == None:
            # If not present, print error.
            print("Could not find \'", name, "\' item.")
        else:
            # Else print out the metaData value.
            print(name, "=\'", metaData, "\'")
    else:
        # Print an error messagge if the file
        # could not be opened.
        print("Could not open the input image file: ", inputFile)

# A function to read a meta-data item
# from a image band

def listImageMetaData(inputFile):
    # Open the dataset in Read Only mode
    dataset = gdal.Open(inputFile, gdal.GA_ReadOnly)
    # Check that the image has been opened.
    if not dataset is None:
        # Get the meta-data dictionary
        metaData = dataset.GetMetaData_Dict()
        # Check it has entries.
        if len(metaData) == 0:
            # if it has no entries return
            # error message.
            print("There is no image meta-data.")
        else:
            # Otherwise, print out the
            # meta-data.
            # Loop through each entry.
            for metaItem in metaData:
                print(metaItem)
    else:
        # Print an error message if the file
        # could not be opened.
        print("Could not open the input image file", inputFile)

# This is the first part of the script to
# be executed.

def copyproj(src_fn, dst_fn, gt=True):
    """Copy projection and geotransform from one raster file to another
    """
    src_ds = gdal.Open(src_fn, gdal.GA_ReadOnly)
    dst_ds = gdal.Open(dst_fn, gdal.GA_Update)
    dst_ds.SetProjection(src_ds.GetProjection())
    if gt:
        src_gt = np.array(src_ds.GetGeoTransform())
        src_dim = np.array([src_ds.RasterXSize, src_ds.RasterYSize])
        dst_dim = np.array([dst_ds.RasterXSize, dst_ds.RasterYSize])
        #This preserves dst_fn resolution
        if np.any(src_dim != dst_dim):
            res_factor = src_dim/dst_dim.astype(float)
            src_gt[[1, 5]] *= max(res_factor)
            #src_gt[[1, 5]] *= min(res_factor)
            #src_gt[[1, 5]] *= res_factor
        dst_ds.SetGeoTransform(src_gt)
    src_ds = None
    dst_ds = None

def rasterize(self,inRaster,inShape,inField):
        filename = tempfile.mktemp('.tif')        
        data = gdal.Open(inRaster,gdal.GA_ReadOnly)
        shp = ogr.Open(inShape)

        lyr = shp.GetLayer()

        driver = gdal.GetDriverByName('GTiff')
        dst_ds = driver.Create(filename,data.RasterXSize,data.RasterYSize, 1,gdal.GDT_Byte)
        dst_ds.SetGeoTransform(data.GetGeoTransform())
        dst_ds.SetProjection(data.GetProjection())
        OPTIONS = 'ATTRIBUTE='+inField
        gdal.RasterizeLayer(dst_ds, [1], lyr, None,options=[OPTIONS])
        data,dst_ds,shp,lyr=None,None,None,None


        return filename

def get_mask(mask_file):
    mask_ds = gdal.Open(mask_file, gdal.GA_ReadOnly)
    mask_data = mask_ds.GetRasterBand(1).ReadAsArray()
    mask = mask_data != MASK_VALUES_TO_KEEP
    mask_ds = None  # close dataset
    return mask

def average(input_dir, out_dir, size):
    input_dir = abspath(input_dir)
    log.info('Reading tifs from {}'.format(input_dir))
    tifs = glob.glob(join(input_dir, '*.tif'))

    for tif in tifs:
        data_source = gdal.Open(tif, gdal.GA_ReadOnly)
        band = data_source.GetRasterBand(1)
        # data_type = gdal.GetDataTypeName(band.DataType)
        data = band.ReadAsArray()
        no_data_val = band.GetNoDataValue()
        averaged_data = filter_data(data, size, no_data_val)
        log.info('Calculated average for {}'.format(basename(tif)))

        output_file = join(out_dir, 'average_' + basename(tif))
        out_ds = gdal.GetDriverByName('GTiff').Create(
            output_file, data_source.RasterXSize, data_source.RasterYSize,
            1, band.DataType)
        out_band = out_ds.GetRasterBand(1)
        out_band.WriteArray(averaged_data)
        out_ds.SetGeoTransform(data_source.GetGeoTransform())
        out_ds.SetProjection(data_source.GetProjection())
        out_band.FlushCache()  # Write data to disc
        out_ds = None  # close out_ds
        data_source = None  # close dataset

        log.info('Finished converting {}'.format(basename(tif)))

def gdalaverage(input_dir, out_dir, size):
    """
    average data using gdal's averaging method.
    Parameters
    ----------
    input_dir: str
        input dir name of the tifs that needs to be averaged
    out_dir: str
        output dir name
    size: int, optional
        size of kernel
    Returns
    -------

    """
    input_dir = abspath(input_dir)
    log.info('Reading tifs from {}'.format(input_dir))
    tifs = glob.glob(join(input_dir, '*.tif'))

    process_tifs = np.array_split(tifs, mpiops.chunks)[mpiops.chunk_index]

    for tif in process_tifs:
        data_set = gdal.Open(tif, gdal.GA_ReadOnly)
        # band = data_set.GetRasterBand(1)
        # data_type = gdal.GetDataTypeName(band.DataType)
        # data = band.ReadAsArray()
        # no_data_val = band.GetNoDataValue()
        # averaged_data = filter_data(data, size, no_data_val)
        log.info('Calculated average for {}'.format(basename(tif)))

        output_file = join(out_dir, 'average_' + basename(tif))
        src_gt = data_set.GetGeoTransform()
        tmp_file = '/tmp/tmp_{}.tif'.format(mpiops.chunk_index)
        resample_cmd = [TRANSLATE] + [tif, tmp_file] + \
            ['-tr', str(src_gt[1]*size), str(src_gt[1]*size)] + \
            ['-r', 'bilinear']
        check_call(resample_cmd)
        rollback_cmd = [TRANSLATE] + [tmp_file, output_file] + \
            ['-tr', str(src_gt[1]), str(src_gt[1])]
        check_call(rollback_cmd)
        log.info('Finished converting {}'.format(basename(tif)))

def get_stats(tif, partitions=100):
    ds = gdal.Open(tif, gdal.GA_ReadOnly)
    number_of_bands = ds.RasterCount

    if ds.RasterCount > 1:
        d = {}
        log.info('Found multibanded geotif {}'.format(basename(tif)))
        for b in range(number_of_bands):
            d['{tif}_{b}'.format(tif=tif, b=b)] = band_stats(ds, tif, b + 1,
                                                             partitions)
        return d
    else:
        log.info('Found single band geotif {}'.format(basename(tif)))
        return {tif: band_stats(ds, tif, 1, partitions)}

def get_numpy_stats(tif, writer):
    ds = gdal.Open(tif, gdal.GA_ReadOnly)
    number_of_bands = ds.RasterCount

    if ds.RasterCount > 1:
        log.info('Found multibanded geotif {}'.format(basename(tif)))
        for b in range(number_of_bands):
            write_rows(stats=number_of_bands(ds, tif, b + 1), writer=writer)
    else:
        log.info('Found single band geotif {}'.format(basename(tif)))
        write_rows(stats=number_of_bands(ds, tif, 1), writer=writer)

def readBandMetaData(inputFile, band, name):
    # Open the dataset in Read Only mode
    dataset = gdal.Open(inputFile, gdal.GA_ReadOnly)
    # Check that the image has opened
    if not dataset is None:
        # Get the image band
        imgBand = dataset.GetRasterBand(band)
        # Check the image band was available
        if not imgBand is None:
            # Get the meta data value specified
            metaData = imgBand.GetMetaDataItem(name)
            # Check that it is present
            if metaData == None:
                # If not present, print error.
                print("Could not find \'", name, "\' item.")
            else:
                # Else print out the metadata value
                print(name, "=\'", metaData, "\'")

        else:
            # Print out an error message.
            print("Could not open the image band: ", band)

    else:
        # Print an error message if the file
        # could not be opened.
        print("Could not open the input image file: ", inputFile)

# A function to read a meta-data item
# from an image

def listBandMetaData(inputFile, band):
    # Open the dataset in Read Only Mode
    dataset = gdal.Open(inputFile, gdal.GA_ReadOnly)
    # Check that the image has been opened.
    if not dataset is None:
        # Get the image band
        imgBand = dataset.GetRasterBand(band)
        # Check the image band was available.
        if not imgBand is None:
            # Get the meta-data dictionary
            metaData = imgBand.GetMetaData_Dict()
            # Check it has entries.
            if len(metaData) == 0:
                # If it has no entries return
                # error message.
                print("There is no image meta-data.")
            else:
                # Otherwise, print out the
                # meta-data.
                # Loop through each entry.
                for metaItem in metaData:
                    print(metaItem)
        else:
            # Print out an error message
            print("Could not open the image bands: ", band)

    else:
        # Print an error message if the file
        # could not be opened.
        print("Could not open the input image file", inputFile)

def fn_getds(fn):
    """Wrapper around gdal.Open()
    """
    ds = None
    if fn_check(fn):
        ds = gdal.Open(fn, gdal.GA_ReadOnly)
    else:
        print("Unable to find %s" % fn)
    return ds

def get_ndv_fn(fn):
    ds = gdal.Open(fn, gdal.GA_ReadOnly)
    return get_ndv_ds(ds)

#Want to modify to handle multi-band images and return list of ndv

def memwarp_multi_fn(src_fn_list, res='first', extent='intersection', t_srs='first', r='cubic', verbose=True, dst_ndv=0):
    """Helper function for memwarp of multiple input filenames
    """
    #Should implement proper error handling here
    if not iolib.fn_list_check(src_fn_list):
        sys.exit('Missing input file(s)')
    src_ds_list = [gdal.Open(fn, gdal.GA_ReadOnly) for fn in src_fn_list]
    return memwarp_multi(src_ds_list, res, extent, t_srs, r, verbose=verbose, dst_ndv=dst_ndv)

def diskwarp_multi_fn(src_fn_list, res='first', extent='intersection', t_srs='first', r='cubic', verbose=True, outdir=None, dst_ndv=None):
    """Helper function for diskwarp of multiple input filenames
    """
    #Should implement proper error handling here
    if not iolib.fn_list_check(src_fn_list):
        sys.exit('Missing input file(s)')
    src_ds_list = [gdal.Open(fn, gdal.GA_ReadOnly) for fn in src_fn_list]
    return diskwarp_multi(src_ds_list, res, extent, t_srs, r, verbose=verbose, outdir=outdir, dst_ndv=dst_ndv)

def dem_geoid(dem_fn):
    out_prefix = os.path.splitext(dem_fn)[0]
    adj_fn = out_prefix +'-adj.tif'
    if not os.path.exists(adj_fn):
        import subprocess
        cmd_args = ["-o", out_prefix, dem_fn]
        cmd = ['dem_geoid'] + cmd_args
        #cmd = 'dem_geoid -o %s %s' % (out_prefix, dem_fn)
        print(' '.join(cmd))
        subprocess.call(cmd, shell=False)
    adj_ds = gdal.Open(adj_fn, gdal.GA_ReadOnly)
    #from pygeotools.lib import iolib
    #return iolib.ds_getma(adj_ds, 1)
    return adj_ds

def _change_segment(self, segment=0, mode=gdal.GA_ReadOnly):
      if segment != self._segment:
        self._dataset = gdal.Open(self.object.GetSubDatasets()[segment][0], mode)
        self._segment = segment

def load(self, mode=gdal.GA_ReadOnly, *args, **kwargs):
      self.object = gdal.Open(self.filename, mode, *args, **kwargs)
      if self.object is None:
        raise Exception('Gdal can not determine driver')
      self._dataset = self.object

def gdaldem_wrapper(fn, product='hs', returnma=True, verbose=True):
    """Wrapper for gdaldem functions

    Note: gdaldem is directly avaialable through API as of GDAL v2.1

    https://trac.osgeo.org/gdal/wiki/rfc59.1_utilities_as_a_library

    This function is no longer necessry, and will eventually be removed.
    """
    #These gdaldem functions should be able to ingest masked array
    #Just write out temporary file, or maybe mem vrt?
    valid_opt = ['hillshade', 'hs', 'slope', 'aspect', 'color-relief', 'TRI', 'TPI', 'roughness']
    try:
        open(fn)
    except IOError:
        print("Unable to open %s" %fn)

    if product not in valid_opt: 
        print("Invalid gdaldem option specified")

    import subprocess
    from pygeotools.lib import iolib
    bma = None
    opts = []
    if product == 'hs' or product == 'hillshade':
        product = 'hillshade'
        #opts = ['-compute_edges',]
        out_fn = os.path.splitext(fn)[0]+'_hs_az315.tif'
    else:
        out_fn = os.path.splitext(fn)[0]+'_%s.tif' % product
    if not os.path.exists(out_fn):
        cmd = ['gdaldem', product]
        cmd.extend(opts)
        cmd.extend(iolib.gdal_opt_co)
        cmd.extend([fn, out_fn])
        if verbose:
            print(' '.join(cmd))
            cmd_opt = {}
        else:
            fnull = open(os.devnull, 'w')
            cmd_opt = {'stdout':fnull, 'stderr':subprocess.STDOUT}
        subprocess.call(cmd, shell=False, **cmd_opt)

    if returnma:
        ds = gdal.Open(out_fn, gdal.GA_ReadOnly)
        bma = iolib.ds_getma(ds, 1)
        return bma 
    else:
        return out_fn

def open_data(filename):
    '''
    The function open and load the image given its name. 
    The type of the data is checked from the file and the scipy array is initialized accordingly.
    Input:
        filename: the name of the file
    Output:
        im: the data cube
        GeoTransform: the geotransform information 
        Projection: the projection information
    '''
    data = gdal.Open(filename,gdal.GA_ReadOnly)
    if data is None:
        print 'Impossible to open '+filename
        exit()
    nc = data.RasterXSize
    nl = data.RasterYSize
    d  = data.RasterCount

    # Get the type of the data
    gdal_dt = data.GetRasterBand(1).DataType
    if gdal_dt == gdal.GDT_Byte:
        dt = 'uint8'
    elif gdal_dt == gdal.GDT_Int16:
        dt = 'int16'
    elif gdal_dt == gdal.GDT_UInt16:
        dt = 'uint16'
    elif gdal_dt == gdal.GDT_Int32:
        dt = 'int32'
    elif gdal_dt == gdal.GDT_UInt32:
        dt = 'uint32'

    elif gdal_dt == gdal.GDT_Float32:
        dt = 'float32'
    elif gdal_dt == gdal.GDT_Float64:
        dt = 'float64'
    elif gdal_dt == gdal.GDT_CInt16 or gdal_dt == gdal.GDT_CInt32 or gdal_dt == gdal.GDT_CFloat32 or gdal_dt == gdal.GDT_CFloat64 :
        dt = 'complex64'
    else:
        print 'Data type unkown'
        exit()

    # Initialize the array
    if d == 1:
        im = sp.empty((nl,nc),dtype=dt) 
    else:
        im = sp.empty((nl,nc,d),dtype=dt) 

    if d == 1:
        im[:,:]=data.GetRasterBand(1).ReadAsArray()
    else :
        for i in range(d):
            im[:,:,i]=data.GetRasterBand(i+1).ReadAsArray()

    GeoTransform = data.GetGeoTransform()
    Projection = data.GetProjection()
    data = None
    return im,GeoTransform,Projection