Python ogr.Open方法代码示例

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def strip_bands(in_raster_path, out_raster_path, bands_to_strip):
    in_raster = gdal.Open(in_raster_path)
    out_raster_band_count = in_raster.RasterCount-len(bands_to_strip)
    out_raster = create_matching_dataset(in_raster, out_raster_path, bands=out_raster_band_count)
    out_raster_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    in_raster_array = in_raster.GetVirtualMemArray()

    bands_to_copy = [band for band in range(in_raster_array.shape[0]) if band not in bands_to_strip]

    out_raster_array[...] = in_raster_array[bands_to_copy, :,:]

    out_raster_array = None
    in_raster_array = None
    out_raster = None
    in_raster = None

    return out_raster_path 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def flatten_probability_image(prob_image, out_path):
    """
    Takes a probability output from classify_image and flattens it into a single layer containing only the maximum
    value from each pixel.

    Parameters
    ----------
    prob_image
        The path to a probability image.
    out_path
        The place to save the flattened image.

    """
    prob_raster = gdal.Open(prob_image)
    out_raster = create_matching_dataset(prob_raster, out_path, bands=1)
    prob_array = prob_raster.GetVirtualMemArray()
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    out_array[:, :] = prob_array.max(axis=0)
    out_array = None
    prob_array = None
    out_raster = None
    prob_raster = None 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def clip_raster_to_intersection(raster_to_clip_path, extent_raster_path, out_raster_path, is_landsat=False):
    """
    Clips one raster to the extent proivded by the other raster, and saves the result at out_raster_path.
    Assumes both raster_to_clip and extent_raster are in the same projection.
    Parameters
    ----------
    raster_to_clip_path
        The location of the raster to be clipped.
    extent_raster_path
        The location of the raster that will provide the extent to clip to
    out_raster_path
        A location for the finished raster
    """

    with TemporaryDirectory() as td:
        temp_aoi_path = os.path.join(td, "temp_clip.shp")
        get_extent_as_shp(extent_raster_path, temp_aoi_path)
        ext_ras = gdal.Open(extent_raster_path)
        proj = osr.SpatialReference(wkt=ext_ras.GetProjection())
        srs_id = int(proj.GetAttrValue('AUTHORITY', 1))
        clip_raster(raster_to_clip_path, temp_aoi_path, out_raster_path, srs_id, flip_x_y = is_landsat) 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def raster_to_array(rst_pth):
    """Reads in a raster file and returns a N-dimensional array.

    Parameters
    ----------
    rst_pth
        Path to input raster.
    Returns
    -------
    As N-dimensional array.
    """
    log = logging.getLogger(__name__)
    in_ds = gdal.Open(rst_pth)
    out_array = in_ds.ReadAsArray()

    return out_array 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def calc_ndvi(raster_path, output_path):
    raster = gdal.Open(raster_path)
    out_raster = create_matching_dataset(raster, output_path, datatype=gdal.GDT_Float32)
    array = raster.GetVirtualMemArray()
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    R = array[2, ...]
    I = array[3, ...]
    out_array[...] = (R-I)/(R+I)

    out_array[...] = np.where(out_array == -2147483648, 0, out_array)

    R = None
    I = None
    array = None
    out_array = None
    raster = None
    out_raster = None 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def apply_image_function(in_paths, out_path, function, out_datatype = gdal.GDT_Int32):
    """Applies a pixel-wise function across every image. Assumes each image is exactly contiguous and, for now,
    single-banded. function() should take a list of values and return a single value."""
    rasters = [gdal.Open(in_path) for in_path in in_paths]
    raster_arrays = [raster.GetVirtualMemArray() for raster in rasters]
    in_array = np.stack(raster_arrays, axis=0)

    out_raster = create_matching_dataset(rasters[0], out_path=out_path, datatype=out_datatype)
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    out_array[...] = np.apply_along_axis(function, 0, in_array)

    # Deallocating. Not taking any chances here.
    out_array = None
    out_raster = None
    in_array = None
    for raster_array, raster in zip(raster_arrays, rasters):
        raster_array = None
        raster = None 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def create_mask_from_model(image_path, model_path, model_clear=0, num_chunks=10, buffer_size=0):
    """Returns a multiplicative mask (0 for cloud, shadow or haze, 1 for clear) built from the model at model_path."""
    from pyeo.classification import classify_image  # Deferred import to deal with circular reference
    with TemporaryDirectory() as td:
        log = logging.getLogger(__name__)
        log.info("Building cloud mask for {} with model {}".format(image_path, model_path))
        temp_mask_path = os.path.join(td, "cat_mask.tif")
        classify_image(image_path, model_path, temp_mask_path, num_chunks=num_chunks)
        temp_mask = gdal.Open(temp_mask_path, gdal.GA_Update)
        temp_mask_array = temp_mask.GetVirtualMemArray()
        mask_path = get_mask_path(image_path)
        mask = create_matching_dataset(temp_mask, mask_path, datatype=gdal.GDT_Byte)
        mask_array = mask.GetVirtualMemArray(eAccess=gdal.GF_Write)
        mask_array[:, :] = np.where(temp_mask_array != model_clear, 0, 1)
        temp_mask_array = None
        mask_array = None
        temp_mask = None
        mask = None
        if buffer_size:
            buffer_mask_in_place(mask_path, buffer_size)
        log.info("Cloud mask for {} saved in {}".format(image_path, mask_path))
        return mask_path 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def create_mask_from_class_map(class_map_path, out_path, classes_of_interest, buffer_size=0, out_resolution=None):
    """Creates a mask from a classification mask: 1 for each pixel containing one of classes_of_interest, otherwise 0"""
    # TODO: pull this out of the above function
    class_image = gdal.Open(class_map_path)
    class_array = class_image.GetVirtualMemArray()
    mask_array = np.isin(class_array, classes_of_interest)
    out_mask = create_matching_dataset(class_image, out_path, datatype=gdal.GDT_Byte)
    out_array = out_mask.GetVirtualMemArray(eAccess=gdal.GA_Update)
    np.copyto(out_array, mask_array)
    class_array = None
    class_image = None
    out_array = None
    out_mask = None
    if out_resolution:
        resample_image_in_place(out_path, out_resolution)
    if buffer_size:
        buffer_mask_in_place(out_path, buffer_size)
    return out_path 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def create_mask_from_fmask(in_l1_dir, out_path):
    log = logging.getLogger(__name__)
    log.info("Creating fmask for {}".format(in_l1_dir))
    with TemporaryDirectory() as td:
        temp_fmask_path = os.path.join(td, "fmask.tif")
        apply_fmask(in_l1_dir, temp_fmask_path)
        fmask_image = gdal.Open(temp_fmask_path)
        fmask_array = fmask_image.GetVirtualMemArray()
        out_image = create_matching_dataset(fmask_image, out_path, datatype=gdal.GDT_Byte)
        out_array = out_image.GetVirtualMemArray(eAccess=gdal.GA_Update)
        log.info("fmask created, converting to binary cloud/shadow mask")
        out_array[:,:] = np.isin(fmask_array, (2, 3, 4), invert=True)
        out_array = None
        out_image = None
        fmask_array = None
        fmask_image = None
        resample_image_in_place(out_path, 10) 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def importAr(to_cur, filename):
    print("Importing:", filename)
    dataSource = ogr.Open(filename)
    dataLayer = dataSource.GetLayer(0)

    for feature in dataLayer:
        geom = feature.GetGeometryRef().ExportToWkt()
        id = feature.GetField("poly_id")
        objtype = feature.GetField("objtype")
        artype = feature.GetField("artype")
        arskogbon = feature.GetField("arskogbon")
        artreslag = feature.GetField("artreslag")
        argrunnf = feature.GetField("argrunnf")

        to_tuple = ( id, objtype, artype, arskogbon, artreslag, argrunnf, geom)
        to_cur.execute("""INSERT into ar_bygg (id, objtype, artype, arskogbon, artreslag, argrunnf, geom)
                        SELECT %s, %s, %s, %s, %s, %s, ST_GeometryFromText(%s);""",
                   to_tuple)

    to_conn.commit()
    dataSource.Destroy() 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def calc_raster_terrain_fixed_elevation(crs, elevation, grid_size, raster_path, locator, x_max, x_min, y_max,
                                        y_min):
    # local variables:
    temp_shapefile = locator.get_temporary_file("terrain.shp")
    cols = int((x_max - x_min) / grid_size)
    rows = int((y_max - y_min) / grid_size)
    shapes = Polygon([[x_min, y_min], [x_max, y_min], [x_max, y_max], [x_min, y_max], [x_min, y_min]])
    geodataframe = Gdf(index=[0], crs=crs, geometry=[shapes])
    geodataframe.to_file(temp_shapefile)
    # 1) opening the shapefile
    source_ds = ogr.Open(temp_shapefile)
    source_layer = source_ds.GetLayer()
    target_ds = gdal.GetDriverByName('GTiff').Create(raster_path, cols, rows, 1, gdal.GDT_Float32)  ##COMMENT 2
    target_ds.SetGeoTransform((x_min, grid_size, 0, y_max, 0, -grid_size))  ##COMMENT 3
    # 5) Adding a spatial reference ##COMMENT 4
    target_dsSRS = osr.SpatialReference()
    target_dsSRS.ImportFromProj4(crs)
    target_ds.SetProjection(target_dsSRS.ExportToWkt())
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(-9999)  ##COMMENT 5
    gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[elevation])  ##COMMENT 6
    target_ds = None  # closing the file 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def add_markers(fmap, json_fn):
    ds = ogr.Open(json_fn)
    lyr = ds.GetLayer()
    for row in lyr:
        geom = row.geometry()
        color = colors[row.GetField('status')]
        fmap.circle_marker([geom.GetY(), geom.GetX()],
                           line_color=color,
                           fill_color=color,
                           radius=5000,
                           popup=get_popup(row.items()))

# Because this version of make_map was defined after the version in
# listing4_3 (since it was already imported and evaluated), then this is
# the version that will be used. Unlike in the book text, the first two
# lines of the function here have been changed to reference get_state_geom,
# save_state_gauges, get_bbox, and get_center from the listing4_3 module. 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def save_state_gauges(out_fn, bbox=None):
    """Save stream gauge data to a geojson file."""
    url = 'http://gis.srh.noaa.gov/arcgis/services/ahps_gauges/' + \
          'MapServer/WFSServer'
    parms = {
        'version': '1.1.0',
        'typeNames': 'ahps_gauges:Observed_River_Stages',
        'srsName': 'urn:ogc:def:crs:EPSG:6.9:4326',
    }
    if bbox:
        parms['bbox'] = bbox
    try:
        request = 'WFS:{0}?{1}'.format(url, urllib.urlencode(parms))
    except:
        request = 'WFS:{0}?{1}'.format(url, urllib.parse.urlencode(parms))
    wfs_ds = ogr.Open(request)
    if wfs_ds is None:
        raise RuntimeError('Could not open WFS.')
    wfs_lyr = wfs_ds.GetLayer(0)

    driver = ogr.GetDriverByName('GeoJSON')
    if os.path.exists(out_fn):
        driver.DeleteDataSource(out_fn)
    json_ds = driver.CreateDataSource(out_fn)
    json_ds.CopyLayer(wfs_lyr, '') 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def open_gpx_file(gpx_file):
    """
    Opens input GPX file
    :param gpx_file: Input GPX file
    :return: File object
    :rtype: Object
    """
    gpx_data_source = ogr.Open(gpx_file)
    if not gpx_data_source:
        raise Exception(
            "File {0} could not be accessed. "
            "It could be in use by another application".format(
                os.path.basename(gpx_file)
            )
        )
    return gpx_data_source 

# 需要导入模块: from osgeo import ogr [as 别名]
# 或者: from osgeo.ogr import Open [as 别名]
def get_feature_layer(gpx_data_source, feature_type):
    """
    Gets valid feature layer from the GPX file based on
    the selected feature type
    :param gpx_data_source: Open file object
    :param feature_type: User feature type input
    :return  gpx_layer: GPX file feature layer
    :return feature_count: Number of features
    :rtype gpx_layer: Layer object
    :rtype feature_count: Integer
    """
    if feature_type >= 0:
        if FEATURE_TYPES[feature_type] == 'track' or \
                        FEATURE_TYPES[feature_type] == 'route':
            feature_type = '{}_points'.format(FEATURE_TYPES[feature_type])
        else:
            feature_type = '{}s'.format(FEATURE_TYPES[feature_type])
        gpx_layer = gpx_data_source.GetLayerByName(feature_type)
        feature_count = gpx_layer.GetFeatureCount()
        return None if feature_count == 0 else gpx_layer, feature_count