Python gdal.GA_ReadOnly方法代码示例

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

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

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

# 需要导入模块: from osgeo import gdal [as 别名]
# 或者: from osgeo.gdal import GA_ReadOnly [as 别名]
def _init_attrs(self, filenames):
        self.filenames = filenames
        self.datasets = [gdal.Open(f, gdal.GA_ReadOnly) for
                         f in self.filenames]

        self.n_image = len(filenames)
        self.n_band = self.datasets[0].RasterCount
        self.n_col = int(self.datasets[0].RasterXSize)
        self.datatype = gdal_array.GDALTypeCodeToNumericTypeCode(
            self.datasets[0].GetRasterBand(1).DataType)

        self.dataset_bands = [
            [ds.GetRasterBand(i + 1) for i in range(self.n_band)]
            for ds in self.datasets
        ] 

# 需要导入模块: from osgeo import gdal [as 别名]
# 或者: from osgeo.gdal import GA_ReadOnly [as 别名]
def loadColisions():
	global obstacle_image, heightMapFile, array, width, height, minHeight, maxHeight, minColisionHeight, maxColisionHeight
	global barrierN, barrierS, barrierE, barrierW, barrierNE, barrierNW, barrierSE, barrierSW, barrier
	# loading height map image
	dataset = gdal.Open(heightMapFile, gdal.GA_ReadOnly)
	try:
		srcband = dataset.GetRasterBand(1)
	except RuntimeError, e:
        	print 'No band %i found' % band_num
        	print e
	        sys.exit(1) 

# 需要导入模块: from osgeo import gdal [as 别名]
# 或者: from osgeo.gdal import GA_ReadOnly [as 别名]
def GetGeoInfo(FileName):
    """This gets information from the raster file using gdal

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

    Return:
        float: A vector that contains:
            * NDV: the nodata values
            * xsize: cellsize in x direction
            * ysize: cellsize in y direction
            * GeoT: the tranform (a string)
            * Projection: the Projection (a string)
            * DataType: The type of data (an int explaing the bits of each data element)

    Author: SMM
    """


    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
#==============================================================================

#==============================================================================
# This gets the UTM zone, if it exists 

# 需要导入模块: from osgeo import gdal [as 别名]
# 或者: from osgeo.gdal import GA_ReadOnly [as 别名]
def PlotRaster(RasterFile, Map, alpha=1.):

    """
    Description goes here...

    MDH
    """

    print "Plotting raster..."

    #Read data
    gdata = gdal.Open(RasterFile, gdal.GA_ReadOnly)
    geo = gdata.GetGeoTransform()
    Data = gdata.ReadAsArray()

    # make topodat a masked array, masking values lower than sea level.
    Data = ma.masked_where(Data < 0, Data)

    #setup meshgrid for raster plotting
    xres = geo[1]
    yres = geo[5]
    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

    x,y = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres]
    x,y = Map(x,y)
    Map.pcolormesh(x, y, Data.T, cmap=plt.cm.Greys, alpha=alpha) 

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

# 需要导入模块: from osgeo import gdal [as 别名]
# 或者: from osgeo.gdal import GA_ReadOnly [as 别名]
def _dataset_from_dir(folder):
        files = set(f for f in os.scandir(folder) if f.is_file())
        for f in files:
            if op.splitext(f.name)[-1] == '':
                break
        else:
            raise ImportError('could not load dataset from %s' % folder)

        return gdal.Open(f.path, gdal.GA_ReadOnly) 

# 需要导入模块: from osgeo import gdal [as 别名]
# 或者: from osgeo.gdal import GA_ReadOnly [as 别名]
def getProjectionSystem(self, filename):
        '''
        Testing with Greenland.
        '''
        if not filename:
            raise Exception('File {0} does not exist'.format(filename))

        from osgeo import gdal, osr
        ds = gdal.Open(filename, gdal.GA_ReadOnly)
        srs = osr.SpatialReference()
        srs.ImportFromWkt(ds.GetProjection())
        srs.AutoIdentifyEPSG()
        ds = None
#        pdb.set_trace()

        if srs.IsGeographic():
            epsgstr = srs.GetAuthorityCode('GEOGCS')
        elif srs.IsProjected():
            epsgstr = srs.GetAuthorityCode('PROJCS')
        elif srs.IsLocal():
            raise Exception('Local coordinate system encountered')
        else:
            raise Exception('Non-standard coordinate system encountered')
        if not epsgstr:  #Empty string->use shell command gdalsrsinfo for last trial
            cmd = 'gdalsrsinfo -o epsg {0}'.format(filename)
            epsgstr = subprocess.check_output(cmd, shell=True)
#            pdb.set_trace()
            epsgstr = re.findall("EPSG:(\d+)", str(epsgstr))[0]
#            pdb.set_trace()
        if not epsgstr:  #Empty string
            raise Exception('Could not auto-identify epsg code')
#        pdb.set_trace()
        epsgcode = int(epsgstr)
#        pdb.set_trace()
        return epsgcode 

# 需要导入模块: from osgeo import gdal [as 别名]
# 或者: from osgeo.gdal import GA_ReadOnly [as 别名]
def getProjectionSystem(self):
        '''
        Testing with Greenland.
        '''
        if not self.demname:
            raise Exception('At least the DEM parameter must be set for geogrid')

        from osgeo import gdal, osr
        ds = gdal.Open(self.demname, gdal.GA_ReadOnly)
        srs = osr.SpatialReference()
        srs.ImportFromWkt(ds.GetProjection())
        srs.AutoIdentifyEPSG()
        ds = None
#        pdb.set_trace()

        if srs.IsGeographic():
            epsgstr = srs.GetAuthorityCode('GEOGCS')
        elif srs.IsProjected():
            epsgstr = srs.GetAuthorityCode('PROJCS')
        elif srs.IsLocal():
            raise Exception('Local coordinate system encountered')
        else:
            raise Exception('Non-standard coordinate system encountered')
        if not epsgstr:  #Empty string->use shell command gdalsrsinfo for last trial
            cmd = 'gdalsrsinfo -o epsg {0}'.format(self.demname)
            epsgstr = subprocess.check_output(cmd, shell=True)
#            pdb.set_trace()
            epsgstr = re.findall("EPSG:(\d+)", str(epsgstr))[0]
#            pdb.set_trace()
        if not epsgstr:  #Empty string
            raise Exception('Could not auto-identify epsg code')
#        pdb.set_trace()
        epsgcode = int(epsgstr)
#        pdb.set_trace()
        return epsgcode 

# 需要导入模块: from osgeo import gdal [as 别名]
# 或者: from osgeo.gdal import GA_ReadOnly [as 别名]
def read_metadata(self, img_name):

        # Read the metadata
        dataset = gdal.Open(img_name, gdal.GA_ReadOnly)
        metadata = dataset.GetMetadata()
        rpc_data = dataset.GetMetadata('RPC')

        # read image and get size
        img = dataset.ReadAsArray()
        self.rows = img.shape[1]
        self.columns = img.shape[2]
        img = None

        # Extract RPC metadata fields
        self.lon_off = float(rpc_data['LONG_OFF'])
        self.lon_scale = float(rpc_data['LONG_SCALE'])
        self.lat_off = float(rpc_data['LAT_OFF'])
        self.lat_scale = float(rpc_data['LAT_SCALE'])
        self.height_off = float(rpc_data['HEIGHT_OFF'])
        self.height_scale = float(rpc_data['HEIGHT_SCALE'])
        self.line_off = float(rpc_data['LINE_OFF'])
        self.line_scale = float(rpc_data['LINE_SCALE'])
        self.samp_off = float(rpc_data['SAMP_OFF'])
        self.samp_scale = float(rpc_data['SAMP_SCALE'])
        self.line_num_coeff = np.asarray(rpc_data['LINE_NUM_COEFF'].split(), dtype=np.float64)
        self.line_den_coeff = np.asarray(rpc_data['LINE_DEN_COEFF'].split(), dtype=np.float64)
        self.samp_num_coeff = np.asarray(rpc_data['SAMP_NUM_COEFF'].split(), dtype=np.float64)
        self.samp_den_coeff = np.asarray(rpc_data['SAMP_DEN_COEFF'].split(), dtype=np.float64)

        # Close the dataset
        dataset = None

    # Forward project normalized WGS84 array to image coordinates
    # Rows and columns are computed with separate calls and coefficient arrays
    # Not currently checking the denominator for zero