Python rasterio.mask方法代码示例

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import mask [as 别名]
def polygonize_raster(infile, outfile, mask_value=1, driver='ESRI Shapefile'):

    with rasterio.open(infile) as src:

        image = src.read(1)

        if mask_value is not None:
            mask = image == mask_value
        else:
            mask = None

        results = (
            {'properties': {'raster_val': v}, 'geometry': s}
            for i, (s, v)
            in enumerate(
                shapes(image, mask=mask, transform=src.transform)))

        with fiona.open(
            outfile, 'w',
            driver=driver,
            crs=src.crs,
            schema={'properties': [('raster_val', 'int')],
                    'geometry': 'Polygon'}) as dst:

            dst.writerecords(results) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import mask [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 mask [as 别名]
def mask_invalid_depth(depth):
    # depth = depth # bigger than the minimum 
    # mask numbers > depth_max
    depth_mask = depth < 9000 # everest has 8.848m
    return depth * depth_mask # values to ignore will be 0 # thats ugly 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import mask [as 别名]
def clip_raster(raster, boundary, boundary_layer=None):
    """Clip the raster to the given administrative boundary.

    Parameters
    ----------
    raster : string, pathlib.Path or rasterio.io.DataSetReader
        Location of or already opened raster.
    boundary : string, pathlib.Path or geopandas.GeoDataFrame
        The polygon by which to clip the raster.
    boundary_layer : string, optional
        For multi-layer files (like GeoPackage), specify the layer to be used.


    Returns
    -------
    tuple
        Three elements:
            clipped : numpy.ndarray
                Contents of clipped raster.
            affine : affine.Affine()
                Information for mapping pixel coordinates
                to a coordinate system.
            crs : dict
                Dict of the form {'init': 'epsg:4326'} defining the coordinate
                reference system of the raster.

    """

    if isinstance(raster, Path):
        raster = str(raster)
    if isinstance(raster, str):
        raster = rasterio.open(raster)

    if isinstance(boundary, Path):
        boundary = str(boundary)
    if isinstance(boundary, str):
        if ".gpkg" in boundary:
            driver = "GPKG"
        else:
            driver = None  # default to shapefile
            boundary_layer = ""  # because shapefiles have no layers

        boundary = gpd.read_file(boundary, layer=boundary_layer, driver=driver)

    if not (boundary.crs == raster.crs or boundary.crs == raster.crs.data):
        boundary = boundary.to_crs(crs=raster.crs)
    coords = [json.loads(boundary.to_json())["features"][0]["geometry"]]

    # mask/clip the raster using rasterio.mask
    clipped, affine = mask(dataset=raster, shapes=coords, crop=True)

    if len(clipped.shape) >= 3:
        clipped = clipped[0]

    return clipped, affine, raster.crs 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import mask [as 别名]
def _predfun(self, img, estimator):
        """Prediction function for classification or regression response.

        Parameters
        ----
        img : tuple (window, numpy.ndarray)
            A window object, and a 3d ndarray of raster data with the dimensions in
            order of (band, rows, columns).

        estimator : estimator object implementing 'fit'
            The object to use to fit the data.

        Returns
        -------
        numpy.ndarray
            2d numpy array representing a single band raster containing the
            classification or regression result.
        """
        window, img = img

        if not isinstance(img, pd.DataFrame):
            # reshape each image block matrix into a 2D matrix
            # first reorder into rows, cols, bands(transpose)
            # then resample into 2D array (rows=sample_n, cols=band_values)
            n_features, rows, cols = img.shape[0], img.shape[1], img.shape[2]
            n_samples = rows * cols
            flat_pixels = img.transpose(1, 2, 0).reshape((n_samples, n_features))

            # create mask for NaN values and replace with number
            flat_pixels_mask = flat_pixels.mask.copy()
            flat_pixels = flat_pixels.filled(0)
        
        else:
            flat_pixels = img
            flat_pixels_mask = pd.isna(flat_pixels).values
            flat_pixels = flat_pixels.fillna(0)
            flat_pixels = flat_pixels.values

        # predict and replace mask
        result_cla = estimator.predict(flat_pixels)
        result_cla = np.ma.masked_array(
            data=result_cla, mask=flat_pixels_mask.any(axis=1)
        )

        # reshape the prediction from a 1D into 3D array [band, row, col]
        result_cla = result_cla.reshape((1, window.height, window.width))

        return result_cla 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import mask [as 别名]
def _probfun(img, estimator):
        """Class probabilities function.

        Parameters
        ----------
        img : tuple (window, numpy.ndarray)
            A window object, and a 3d ndarray of raster data with the dimensions in
            order of (band, rows, columns).

        estimator : estimator object implementing 'fit'
            The object to use to fit the data.

        Returns
        -------
        numpy.ndarray
            Multi band raster as a 3d numpy array containing the probabilities
            associated with each class. ndarray dimensions are in the order of
            (class, row, column).
        """
        window, img = img

        if not isinstance(img, pd.DataFrame):
            # reshape each image block matrix into a 2D matrix
            # first reorder into rows, cols, bands (transpose)
            # then resample into 2D array (rows=sample_n, cols=band_values)
            n_features, rows, cols = img.shape[0], img.shape[1], img.shape[2]
            mask2d = img.mask.any(axis=0)
            n_samples = rows * cols
            flat_pixels = img.transpose(1, 2, 0).reshape((n_samples, n_features))
            flat_pixels = flat_pixels.filled(0)
        
        else:
            flat_pixels = img
            mask2d = pd.isna(flat_pixels).values
            mask2d = mask2d.reshape((window.height, window.width, flat_pixels.shape[1]))
            mask2d = mask2d.any(axis=2)
            flat_pixels = flat_pixels.fillna(0)
            flat_pixels = flat_pixels.values

        # predict probabilities
        result_proba = estimator.predict_proba(flat_pixels)

        # reshape class probabilities back to 3D image [iclass, rows, cols]
        result_proba = result_proba.reshape((window.height, window.width, result_proba.shape[1]))

        # reshape band into rasterio format [band, row, col]
        result_proba = result_proba.transpose(2, 0, 1)

        # repeat mask for n_bands
        mask3d = np.repeat(
            a=mask2d[np.newaxis, :, :], repeats=result_proba.shape[0], axis=0
        )

        # convert proba to masked array
        result_proba = np.ma.masked_array(result_proba, mask=mask3d, fill_value=np.nan)

        return result_proba 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import mask [as 别名]
def _predfun_multioutput(img, estimator):
        """Multi-target prediction function.

        Parameters
        ----------
        img : tuple (window, numpy.ndarray)
            A window object, and a 3d ndarray of raster data with the dimensions in
            order of (band, rows, columns).

        estimator : estimator object implementing 'fit'
            The object to use to fit the data.

        Returns
        -------
        numpy.ndarray
            3d numpy array representing the multi-target prediction result with the
            dimensions in the order of (target, row, column).
        """
        window, img = img

        if not isinstance(img, pd.DataFrame):
            # reshape each image block matrix into a 2D matrix
            # first reorder into rows, cols, bands(transpose)
            # then resample into 2D array (rows=sample_n, cols=band_values)
            n_features, rows, cols = img.shape[0], img.shape[1], img.shape[2]
            mask2d = img.mask.any(axis=0)
            n_samples = rows * cols
            flat_pixels = img.transpose(1, 2, 0).reshape((n_samples, n_features))
            flat_pixels = flat_pixels.filled(0)
        
        else:
            flat_pixels = img
            mask2d = pd.isna(flat_pixels).values
            mask2d = mask2d.reshape((window.height, window.width, flat_pixels.shape[1]))
            mask2d = mask2d.any(axis=2)
            flat_pixels = flat_pixels.fillna(0)
            flat_pixels = flat_pixels.values

        # predict probabilities
        result = estimator.predict(flat_pixels)

        # reshape class probabilities back to 3D image [iclass, rows, cols]
        result = result.reshape((window.height, window.width, result.shape[1]))

        # reshape band into rasterio format [band, row, col]
        result = result.transpose(2, 0, 1)

        # repeat mask for n_bands
        mask3d = np.repeat(a=mask2d[np.newaxis, :, :], repeats=result.shape[0], axis=0)

        # convert proba to masked array
        result = np.ma.masked_array(result, mask=mask3d, fill_value=np.nan)

        return result 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import mask [as 别名]
def load_rgb_and_label(input_list, phase, target_path=None, dfc_preprocessing=0):
    """
    load all depths/labels to label_cache with crop
    """
    # ToDo: fix this mess of having two almost similar functions... 
    import rasterio
    from rasterio.mask import mask
    from rasterio.plot import reshape_as_image
    import geopandas as gpd
    from shapely.geometry import box
    from scipy.misc import imresize
    
    if phase != 'test':
        raster_target = [rasterio.open(path) for path in target_path]
    
    print('Loading {} patches...'.format(phase))
    for i, img_input_path in enumerate(tqdm(input_list)):
        with rasterio.open(img_input_path) as raster_rgb:

            # normalize data between -1 and 1 and append to cache
            raster_rgb_numpy = raster_rgb.read()
            raster_rgb_norm = ((raster_rgb_numpy.astype('float32') / 255.0) * 2.0) - 1.0
            pil_shape = (raster_rgb_norm.shape[-2:])[::-1]
            
            if dfc_preprocessing == 2: # make input same resolution as output
                pil_shape = [dim // 10 for dim in pil_shape] # Resize
                raster_rgb_norm = resize_rgb(raster_rgb_norm, pil_shape, Image.BILINEAR)

            if phase == 'test':
                yield raster_rgb_norm #, pil_depth_patch_shape # To add labels for test
            else:
                labels_patches = []

                # get bounds and crop from depth_image_raster
                bounds = raster_rgb.bounds
                bbox = box(bounds.left, bounds.bottom, bounds.right, bounds.top)
                geo = gpd.GeoDataFrame({'geometry': bbox}, index=[0])
                coords = getFeatures(geo)
                labels_patches_masks = [mask(image_raster, shapes=coords, crop=True) for image_raster in raster_target]

                sem_label = labels_patches_masks[0][0]
            
                if dfc_preprocessing == 0:
                    sem_label = [resize(sem_label[0], pil_shape, Image.NEAREST)]
                elif dfc_preprocessing == 1:
                    pil_shape = (raster_rgb_norm.shape[-2:])[::-1]
                    pil_shape = [dim // 2 for dim in pil_shape] # Resize
                    raster_rgb_norm = resize_rgb(raster_rgb_norm, pil_shape, Image.BILINEAR)
                
                # labels_patches.append(sem_label)
                yield raster_rgb_norm, sem_label 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import mask [as 别名]
def mask_by_shape(infile, outfile, shapefile, to_db=False, datatype='float32',
                  rescale=True, min_value=0.000001, max_value=1, ndv=None,
                  description=True):

    # import shapefile geometries
    with fiona.open(shapefile, 'r') as file:
        features = [feature['geometry'] for feature in file
                    if feature['geometry']]

    # import raster
    with rasterio.open(infile) as src:
        out_image, out_transform = rasterio.mask.mask(src, features, crop=True)
        out_meta = src.meta.copy()
        out_image = np.ma.masked_where(out_image == ndv, out_image)

    
    # unmask array
    out_image = out_image.data
    
    # if to decibel should be applied
    if to_db is True:
        out_image = convert_to_db(out_image)

    if rescale:
        # if we scale to another d
        if datatype != 'float32':

            if datatype == 'uint8':
                out_image = scale_to_int(out_image, min_value, max_value, 'uint8')
            elif datatype == 'uint16':
                out_image = scale_to_int(out_image, min_value, max_value, 'uint16')

    out_meta.update({'driver': 'GTiff', 'height': out_image.shape[1],
                     'width': out_image.shape[2], 'transform': out_transform,
                     'nodata': ndv, 'dtype': datatype, 'tiled': True,
                     'blockxsize': 128, 'blockysize': 128})

    with rasterio.open(outfile, 'w', **out_meta) as dest:
        dest.write(out_image)
        if description:
            dest.update_tags(1, 
                    BAND_NAME='{}'.format(os.path.basename(infile)[:-4]))
            dest.set_band_description(1, 
                    '{}'.format(os.path.basename(infile)[:-4])) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import mask [as 别名]
def save_tile(self, tile_data, mask, profile, dest_fname_base=None):
        """Save a tile created by ``Tiler.tile_generator()``."""
        if dest_fname_base is None:
            dest_fname_root = os.path.splitext(
                os.path.split(self.src_path)[1])[0]
        else:
            dest_fname_root = dest_fname_base
        if self.proj_unit not in ['meter', 'metre']:
            dest_fname = '{}_{}_{}.tif'.format(
                dest_fname_root,
                np.round(profile['transform'][2], 3),
                np.round(profile['transform'][5], 3))
        else:
            dest_fname = '{}_{}_{}.tif'.format(
                dest_fname_root,
                int(profile['transform'][2]),
                int(profile['transform'][5]))
        # if self.cog_output:
        #     dest_path = os.path.join(self.dest_dir, 'tmp.tif')
        # else:
        dest_path = os.path.join(self.dest_dir, dest_fname)

        with rasterio.open(dest_path, 'w',
                           **profile) as dest:
            if profile['count'] == 1:
                dest.write(tile_data[0, :, :], 1)
            else:
                for band in range(1, profile['count'] + 1):
                    # base-1 vs. base-0 indexing...bleh
                    dest.write(tile_data[band-1, :, :], band)
            if self.alpha:
                # write the mask if there's an alpha band
                dest.write(mask, profile['count'] + 1)

            dest.close()

        return dest_path

        # if self.cog_output:
        #     self._create_cog(os.path.join(self.dest_dir, 'tmp.tif'),
        #                      os.path.join(self.dest_dir, dest_fname))
        #     os.remove(os.path.join(self.dest_dir, 'tmp.tif')) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import mask [as 别名]
def fill_all_nodata(self, nodata_fill):
        """
        Fills all tile nodata values with a fill value.
        
        The standard workflow is to run this function only after generating label masks and using the original output 
        from the raster tiler to filter out label pixels that overlap nodata pixels in a tile. For example, 
        solaris.vector.mask.instance_mask will filter out nodata pixels from a label mask if a reference_im is provided,
        and after this step nodata pixels may be filled by calling this method.
        
        nodata_fill : int, float, or str, optional
            Default is to not fill any nodata values. Otherwise, pixels outside of the aoi_boundary and pixels inside 
            the aoi_boundary with the nodata value will be filled. "mean" will fill pixels with the channel-wise mean. 
            Providing an int or float will fill pixels in all channels with the provided value.
            
        Returns: list
            The fill values, in case the mean of the src image should be used for normalization later.
        """
        src = _check_rasterio_im_load(self.src_name)
        if nodata_fill == "mean":
            arr = src.read()
            arr_nan = np.where(arr!=src.nodata, arr, np.nan)
            fill_values = np.nanmean(arr_nan, axis=tuple(range(1, arr_nan.ndim)))
            print('Fill values set to {}'.format(fill_values))
        elif isinstance(nodata_fill, (float, int)):
            fill_values = src.meta['count'] * [nodata_fill]
            print('Fill values set to {}'.format(fill_values))
        else:
            raise TypeError('nodata_fill must be "mean", int, or float. {} was supplied.'.format(nodata_fill))
        src.close()
        for tile_path in self.tile_paths:
            tile_src = rasterio.open(tile_path, "r+")
            tile_data = tile_src.read()
            for i in np.arange(tile_data.shape[0]):
                tile_data[i,...][tile_data[i,...] == tile_src.nodata] = fill_values[i] # set fill value for each band
            if tile_src.meta['count'] == 1:
                tile_src.write(tile_data[0, :, :], 1)
            else:
                for band in range(1, tile_src.meta['count'] + 1):
                    # base-1 vs. base-0 indexing...bleh
                    tile_src.write(tile_data[band-1, :, :], band)
            tile_src.close()
        return fill_values