Python rasterio.features方法代码示例

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def projectShapes(features, toCRS):
    import pyproj
    from functools import partial
    import fiona.crs as fcrs
    from shapely.geometry import shape, mapping
    from shapely.ops import transform as shpTrans

    project = partial(
        pyproj.transform,
        pyproj.Proj(fcrs.from_epsg(4326)),
        pyproj.Proj(toCRS))

    return list(
        {'geometry': mapping(
            shpTrans(
                project,
                shape(feat['geometry']))
        )} for feat in features
        ) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def __init__(self, raster_mask, no_data_value=None, ignore_labels=None):
        if ignore_labels is None:
            ignore_labels = []

        self.geometries = [{'label': int(label),
                            'polygon': Polygon(LinearRing(shp['coordinates'][0]),
                                               [LinearRing(pts) for pts in shp['coordinates'][1:]])}
                           for index, (shp, label) in enumerate(rasterio.features.shapes(raster_mask, mask=None))
                           if (int(label) is not no_data_value) and (int(label) not in ignore_labels)]

        self.areas = np.asarray([entry['polygon'].area for entry in self.geometries])
        self.decomposition = [shapely.ops.triangulate(entry['polygon']) for entry in self.geometries]

        self.label2cc = collections.defaultdict(list)
        for index, entry in enumerate(self.geometries):
            self.label2cc[entry['label']].append(index) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def rasterize_overlapped(self, shapes, out, **rasterize_args):
        """ Rasterize overlapped classes.

        :param shapes: Shapes to be rasterized.
        :type shapes: an iterable of pairs (rasterio.polygon, int)
        :param out: A numpy array to which to rasterize polygon classes.
        :type out: numpy.ndarray
        :param rasterize_args: Keyword arguments to be passed to `rasterio.features.rasterize`.
        :type rasterize_args: dict
        """
        rasters = [rasterio.features.rasterize([shape], out=np.copy(out), **rasterize_args) for shape in shapes]

        overlap_mask = np.zeros(out.shape, dtype=np.bool)
        no_data = self.no_data_value

        out[:] = rasters[0][:]
        for raster in rasters[1:]:
            overlap_mask[(out != no_data) & (raster != no_data) & (raster != out)] = True
            out[raster != no_data] = raster[raster != no_data]

        out[overlap_mask] = self.overlap_value 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def polygonize(self, data=None, mask=None, connectivity=4, transform=None):
        """
        Yield (polygon, value) for each set of adjacent pixels of the same value.
        Wrapper around rasterio.features.shapes

        From rasterio documentation:

        Parameters
        ----------
        data : numpy ndarray
        mask : numpy ndarray
               Values of False or 0 will be excluded from feature generation.
        connectivity : 4 or 8 (int)
                       Use 4 or 8 pixel connectivity.
        transform : affine.Affine
                    Transformation from pixel coordinates of `image` to the
                    coordinate system of the input `shapes`.
        """
        if not _HAS_RASTERIO:
            raise ImportError('Requires rasterio module')
        if data is None:
            data = self.mask.astype(np.uint8)
        if mask is None:
            mask = self.mask
        if transform is None:
            transform = self.affine
        shapes = rasterio.features.shapes(data, mask=mask, connectivity=connectivity,
                                          transform=transform)
        return shapes 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def clip_rasters(folder_in, folder_out, aoi_in, debug=False):
    """Read continental rasters one at a time, clip to AOI and save

    Parameters
    ----------
    folder_in : str, Path
        Path to directory containing rasters.
    folder_out : str, Path
        Path to directory to save clipped rasters.
    aoi_in : str, Path
        Path to an AOI file (readable by Fiona) to use for clipping.
    """

    if isinstance(aoi_in, gpd.GeoDataFrame):
        aoi = aoi_in
    else:
        aoi = gpd.read_file(aoi_in)

    coords = [json.loads(aoi.to_json())["features"][0]["geometry"]]

    for file_path in os.listdir(folder_in):
        if file_path.endswith(".tif"):
            if debug:
                print(f"Doing {file_path}")
            ntl_rd = rasterio.open(os.path.join(folder_in, file_path))
            ntl, affine = mask(dataset=ntl_rd, shapes=coords, crop=True, nodata=0)

            if ntl.ndim == 3:
                ntl = ntl[0]

            save_raster(folder_out / file_path, ntl, affine) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def merge_rasters(folder, percentile=70):
    """Merge a set of monthly rasters keeping the nth percentile value.

    Used to remove transient features from time-series data.

    Parameters
    ----------
    folder : str, Path
        Folder containing rasters to be merged.
    percentile : int, optional (default 70.)
        Percentile value to use when merging using np.nanpercentile.
        Lower values will result in lower values/brightness.

    Returns
    -------
    raster_merged : numpy array
        The merged array.
    affine : affine.Affine
        The affine transformation for the merged raster.
    """

    affine = None
    rasters = []

    for file in os.listdir(folder):
        if file.endswith(".tif"):
            ntl_rd = rasterio.open(os.path.join(folder, file))
            rasters.append(ntl_rd.read(1))

            if not affine:
                affine = ntl_rd.transform

    raster_arr = np.array(rasters)

    raster_merged = np.percentile(raster_arr, percentile, axis=0)

    return raster_merged, affine 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def preprocess_data(tile_size, dataset, only_cache=False):
    """Create features and labels for a given dataset. The features are tiles which contain
    the three RGB bands of the satellite image, so they have the form (tile_size, tile_size, 3).
    Labels are bitmaps with 1 indicating that the corresponding pixel in the satellite image
    represents water."""

    print('_' * 100)
    print("Start preprocessing data.")

    features_train, labels_train = extract_features_and_labels(
        dataset["train"], tile_size, only_cache)
    features_test, labels_test = extract_features_and_labels(
        dataset["test"], tile_size, only_cache)

    return features_train, features_test, labels_train, labels_test 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def extract_features_and_labels(dataset, tile_size, only_cache=False):
    """For each satellite image and its corresponding shapefiles in the dataset create
    tiled features and labels."""
    features = []
    labels = []

    for geotiff_path, shapefile_paths in dataset:
        tiled_features, tiled_labels = create_tiled_features_and_labels(
            geotiff_path, shapefile_paths, tile_size, only_cache)

        features += tiled_features
        labels += tiled_labels

    return features, labels 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def get_geometry_mask(width: int, height: int,
                      geometry: GeometryLike,
                      x_min: float, y_min: float, res: float) -> np.ndarray:
    geometry = convert_geometry(geometry)
    # noinspection PyTypeChecker
    transform = affine.Affine(res, 0.0, x_min,
                              0.0, -res, y_min + res * height)
    return rasterio.features.geometry_mask([geometry],
                                           out_shape=(height, width),
                                           transform=transform,
                                           all_touched=True,
                                           invert=True) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def mask_to_poly(mask, min_polygon_area_th=MIN_AREA):
    shapes = rasterio.features.shapes(mask.astype(np.int16), mask > 0)
    poly_list = []
    mp = shapely.ops.cascaded_union(
        shapely.geometry.MultiPolygon([
            shapely.geometry.shape(shape)
            for shape, value in shapes
        ]))

    if isinstance(mp, shapely.geometry.Polygon):
        df = pd.DataFrame({
            'area_size': [mp.area],
            'poly': [mp],
        })
    else:
        df = pd.DataFrame({
            'area_size': [p.area for p in mp],
            'poly': [p for p in mp],
        })

    df = df[df.area_size > min_polygon_area_th].sort_values(
        by='area_size', ascending=False)
    df.loc[:, 'wkt'] = df.poly.apply(lambda x: shapely.wkt.dumps(
        x, rounding_precision=0))
    df.loc[:, 'bid'] = list(range(1, len(df) + 1))
    df.loc[:, 'area_ratio'] = df.area_size / df.area_size.max()
    return df 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def getBounds(features):
    xy = np.vstack(list(f['geometry']['coordinates'][0] for f in features))
    return coords.BoundingBox(
        xy[:,0].min(),
        xy[:,1].min(),
        xy[:,0].max(),
        xy[:,1].max()
        ) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def getGJSONinfo(geoJSONinfo):
    """
    Loads a lattice of GeoJSON, bounds, and creates a list mapping an on-the-fly UID w/ the actual index value.
    """
    features = list(i for i in filterBadJSON(geoJSONinfo))
    UIDs = list(feat['properties']['qt'] for feat in features)

    featDimensions = int(np.sqrt(len(features)/2.0))
    
    return features, UIDs, featDimensions 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def fillFacets(geoJSONpath, rasterPath, noProject, output, bands, zooming, batchprint, outputGeom, color):
    geoJSON, uidMap, featDims = getGJSONinfo(geoJSONpath)

    rasCRS, rasBounds, rasBands = getRasterInfo(rasterPath)

    bands = handleBandArgs(bands, rasBands)

    if rasCRS['proj'] == 'longlat' or noProject:
        noProject = True
        bounds = getBounds(geoJSON)
    else:
        ogeoJson = geoJSON
        geoJSON = projectShapes(geoJSON, rasCRS)
        bounds = getBounds(geoJSON)

    rasArr, oaff = loadRaster(rasterPath, bands, bounds)

    if min(rasArr.shape[0:2]) < 2 * featDims or zooming:
        rasArr = upsampleRaster(rasArr, featDims, zooming)


    if noProject:
        sampleVals = getRasterValues(geoJSON, rasArr, uidMap, bounds, outputGeom, bands, color)
    else:
        sampleVals = getRasterValues(geoJSON, rasArr, uidMap, bounds, outputGeom, bands, color, ogeoJson)

    if batchprint and outputGeom != True:
        sampleVals = batchStride(sampleVals, int(batchprint))

    if output:
        with open(output, 'w') as oFile:
            oFile.write(json.dumps({
                "type": "FeatureCollection",
                "features": list(sampleVals)
                }))
    else:
        try:
            for feat in sampleVals:
                click.echo(json.dumps(feat).rstrip())
        except IOError as e:
            pass 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def test_compute_metadata_approximate(nodata_type, big_raster_file_nodata, big_raster_file_mask):
    from terracotta.drivers.raster_base import RasterDriver

    if nodata_type == 'nodata':
        raster_file = big_raster_file_nodata
    elif nodata_type == 'masked':
        raster_file = big_raster_file_mask

    with rasterio.open(str(raster_file)) as src:
        data = src.read(1, masked=True)
        valid_data = data.compressed()
        dataset_shape = list(rasterio.features.dataset_features(
            src, bidx=1, band=False, as_mask=True, geographic=True
        ))

    convex_hull = MultiPolygon([shape(s['geometry']) for s in dataset_shape]).convex_hull

    # compare
    mtd = RasterDriver.compute_metadata(str(raster_file), max_shape=(512, 512))

    np.testing.assert_allclose(mtd['valid_percentage'], 100 * valid_data.size / data.size, atol=1)
    np.testing.assert_allclose(
        mtd['range'], (valid_data.min(), valid_data.max()), atol=valid_data.max() / 100
    )
    np.testing.assert_allclose(mtd['mean'], valid_data.mean(), rtol=0.02)
    np.testing.assert_allclose(mtd['stdev'], valid_data.std(), rtol=0.02)

    np.testing.assert_allclose(
        mtd['percentiles'],
        np.percentile(valid_data, np.arange(1, 100)),
        atol=valid_data.max() / 100, rtol=0.02
    )

    assert geometry_mismatch(shape(mtd['convex_hull']), convex_hull) < 0.05 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import features [as 别名]
def test_compute_metadata_nocrick(big_raster_file_nodata, monkeypatch):
    with rasterio.open(str(big_raster_file_nodata)) as src:
        data = src.read(1, masked=True)
        valid_data = data.compressed()
        dataset_shape = list(rasterio.features.dataset_features(
            src, bidx=1, band=False, as_mask=True, geographic=True
        ))

    convex_hull = MultiPolygon([shape(s['geometry']) for s in dataset_shape]).convex_hull

    from terracotta import exceptions
    import terracotta.drivers.raster_base

    with monkeypatch.context() as m:
        m.setattr(terracotta.drivers.raster_base, 'has_crick', False)

        with pytest.warns(exceptions.PerformanceWarning):
            mtd = terracotta.drivers.raster_base.RasterDriver.compute_metadata(
                str(big_raster_file_nodata), use_chunks=True
            )

    # compare
    np.testing.assert_allclose(mtd['valid_percentage'], 100 * valid_data.size / data.size)
    np.testing.assert_allclose(mtd['range'], (valid_data.min(), valid_data.max()))
    np.testing.assert_allclose(mtd['mean'], valid_data.mean())
    np.testing.assert_allclose(mtd['stdev'], valid_data.std())

    # allow error of 1%, since we only compute approximate quantiles
    np.testing.assert_allclose(
        mtd['percentiles'],
        np.percentile(valid_data, np.arange(1, 100)),
        rtol=2e-2
    )

    assert geometry_mismatch(shape(mtd['convex_hull']), convex_hull) < 1e-8