Python rasterio.transform方法代码示例

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def __init__(
        self, in_tile=None, in_data=None, out_profile=None, out_tile=None, tags=None
    ):
        """Prepare data & profile."""
        out_tile = out_tile or in_tile
        validate_write_window_params(in_tile, out_tile, in_data, out_profile)
        self.data = extract_from_array(
            in_raster=in_data,
            in_affine=in_tile.affine,
            out_tile=out_tile
        )
        # use transform instead of affine
        if "affine" in out_profile:
            out_profile["transform"] = out_profile.pop("affine")
        self.profile = out_profile
        self.tags = tags 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def bounds_to_ranges(out_bounds=None, in_affine=None, in_shape=None):
    """
    Return bounds range values from geolocated input.

    Parameters
    ----------
    out_bounds : tuple
        left, bottom, right, top
    in_affine : Affine
        input geolocation
    in_shape : tuple
        input shape

    Returns
    -------
    minrow, maxrow, mincol, maxcol
    """
    return itertools.chain(
        *from_bounds(
            *out_bounds, transform=in_affine, height=in_shape[-2], width=in_shape[-1]
        ).round_lengths(pixel_precision=0).round_offsets(pixel_precision=0).toranges()
    ) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def _reproject(self, eopatch, src_raster):
        """
        Reprojects the raster data from Geopedia's CRS (POP_WEB) to EOPatch's CRS.
        """
        height, width = src_raster.shape

        dst_raster = np.ones((height, width), dtype=self.raster_dtype)

        src_bbox = eopatch.bbox.transform(CRS.POP_WEB)
        src_transform = rasterio.transform.from_bounds(*src_bbox, width=width, height=height)

        dst_bbox = eopatch.bbox
        dst_transform = rasterio.transform.from_bounds(*dst_bbox, width=width, height=height)

        rasterio.warp.reproject(src_raster, dst_raster,
                                src_transform=src_transform, src_crs={'init': CRS.ogc_string(CRS.POP_WEB)},
                                src_nodata=0,
                                dst_transform=dst_transform, dst_crs={'init': CRS.ogc_string(eopatch.bbox.crs)},
                                dst_nodata=self.no_data_val)

        return dst_raster 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def threshold(dists_in, cutoff=0.0):
    """Convert distance array into binary array of connected locations.

    Parameters
    ----------
    dists_in : path-like or numpy array
        2D array output from gridfinder algorithm.
    cutoff : float, optional (default 0.5.)
        Cutoff value below which consider the cells to be grid.

    Returns
    -------
    guess : numpy array
        Binary representation of input array.
    affine: affine.Affine
        Affine transformation for raster.
    """
    if isinstance(dists_in, (str, Path)):
        dists_rd = rasterio.open(dists_in)
        dists_r = dists_rd.read(1)
        affine = dists_rd.transform

        guess = dists_r.copy()
        guess[dists_r > cutoff] = 0
        guess[dists_r <= cutoff] = 1

        return guess, affine

    elif isinstance(dists_in, np.ndarray):
        guess = dists_in.copy()
        guess[dists_in > cutoff] = 0
        guess[dists_in <= cutoff] = 1

        return guess

    else:
        raise ValueError 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def thin(guess_in):
    """
    Use scikit-image skeletonize to 'thin' the guess raster.

    Parameters
    ----------
    guess_in : path-like or 2D array
        Output from threshold().

    Returns
    -------
    guess_skel : numpy array
        Thinned version.
    affine : Affine
        Only if path-like supplied.
    """

    if isinstance(guess_in, (str, Path)):
        guess_rd = rasterio.open(guess_in)
        guess_arr = guess_rd.read(1)
        affine = guess_rd.transform

        guess_skel = skeletonize(guess_arr)
        guess_skel = guess_skel.astype("int32")

        return guess_skel, affine

    elif isinstance(guess_in, np.ndarray):
        guess_skel = skeletonize(guess_in)
        guess_skel = guess_skel.astype("int32")

        return guess_skel

    else:
        raise ValueError 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def run(self, processes=4):
        """TODO"""
        if processes == 1:
            self.pool = MockTub(init_worker, (self.inpaths, self.global_args))
        else:
            self.pool = Pool(processes, init_worker, (self.inpaths, self.global_args))

        self.options["transform"] = guard_transform(self.options["transform"])

        if self.mode == "manual_read":
            reader_worker = manual_reader(self.run_function)
        elif self.mode == "array_read":
            reader_worker = array_reader(self.run_function)
        else:
            reader_worker = simple_reader(self.run_function)

        if isinstance(self.outpath_or_dataset, rasterio.io.DatasetWriter):
            destination = self.outpath_or_dataset
        else:
            destination = rasterio.open(self.outpath_or_dataset, "w", **self.options)

        # Open an output file, work through the function in parallel,
        # and write out the data.
        with destination as dst:
            for data, window in self.pool.imap_unordered(reader_worker, self.windows):
                dst.write(data, window=window)

        self.pool.close()
        self.pool.join() 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def _check_alignment(layers):
        """Check that a list of raster datasets are aligned with the same pixel
        dimensions and geotransforms.

        Parameters
        ----------
        layers : list
            List of pyspatialml.RasterLayer objects.

        Returns
        -------
        dict or False
            Dict of metadata if all layers are spatially aligned, otherwise
            returns False.
        """

        src_meta = []
        for layer in layers:
            src_meta.append(layer.ds.meta.copy())

        if not all(i["crs"] == src_meta[0]["crs"] for i in src_meta):
            Warning(
                "crs of all rasters does not match, " "possible unintended consequences"
            )

        if not all(
            [
                i["height"] == src_meta[0]["height"]
                or i["width"] == src_meta[0]["width"]
                or i["transform"] == src_meta[0]["transform"]
                for i in src_meta
            ]
        ):
            return False

        else:
            return src_meta[0] 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def save_merged_rasters(self, datatype, fileroot=None):
        import rasterio
        from rasterio.merge import merge
        from rasterio.plot import show
        from os.path import join
        import argparse
        import glob
        
        if fileroot == None:
            fileroot = datatype

        root = '{}/{}/{}*.tif'.format(self.save_samples_path, datatype, fileroot)
        filename = '{}/{}/merged_{}.tif'.format(self.save_samples_path, datatype, fileroot)

        files = glob.glob(join(root))
        mosaic_rasters = [rasterio.open(file) for file in files]

        mosaic, out_transform = merge(mosaic_rasters)

        meta = (rasterio.open(files[0])).meta

        meta.update({"driver": "GTif",
                    "height": mosaic.shape[1],
                    "width":  mosaic.shape[2],
                    "transform": out_transform})

        with rasterio.open(filename, "w", **meta) as dest:
            dest.write(mosaic)

        filename = '{}/{}/{}_merged.png'.format(self.save_samples_path, datatype, datatype) 
        self.save_raster_png(mosaic, filename) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def save_merged_rasters(self, datatype, fileroot=None):
        import rasterio
        from rasterio.merge import merge
        from rasterio.plot import show
        from os.path import join
        import argparse
        import glob
        
        if fileroot == None:
            fileroot = datatype

        root = '{}/{}/{}*.tif'.format(self.save_samples_path, datatype, fileroot)
        filename = '{}/{}/merged_{}.tif'.format(self.save_samples_path, datatype, fileroot)

        files = glob.glob(join(root))
        mosaic_rasters = [rasterio.open(file) for file in files]

        mosaic, out_transform = merge(mosaic_rasters)

        meta = (rasterio.open(files[0])).meta

        meta.update({"driver": "GTiff",
                    "height": mosaic.shape[1],
                    "width":  mosaic.shape[2],
                    "transform": out_transform})

        with rasterio.open(filename, "w", **meta) as dest:
            dest.write(mosaic)

        filename = '{}/{}/{}_merged.png'.format(self.save_samples_path, datatype, datatype) 
        self.save_raster_png(mosaic, filename)

        if 'output' in filename or 'target' in filename:
            self.save_height_colormap(filename, mosaic) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def __init__(self, path, cloud):
        import rasterio

        self.in_raster = rasterio.open(path)

        # Check cell size
        cell_size_x, cell_size_y = (
            self.in_raster.transform[0],
            abs(self.in_raster.transform[4]),
        )
        if cell_size_x != cell_size_y:
            print("Cell sizes not equal of input raster, not supported.")
            raise ValueError
        else:
            cell_size = cell_size_x

        self.cloud = cloud
        self.cell_size = cell_size

        min_x, max_x = self.in_raster.bounds[0], self.in_raster.bounds[2]
        min_y, max_y = self.in_raster.bounds[1], self.in_raster.bounds[3]

        self.m = self.in_raster.height
        self.n = self.in_raster.width

        # Create bins
        bins_x = np.searchsorted(
            np.linspace(min_x, max_x, self.n), self.cloud.data.points["x"]
        )
        bins_y = np.searchsorted(
            np.linspace(min_y, max_y, self.m), self.cloud.data.points["y"]
        )

        self.cloud.data.points["bins_x"] = bins_x
        self.cloud.data.points["bins_y"] = bins_y
        self.cells = self.cloud.data.points.groupby(["bins_x", "bins_y"]) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def pansharpen(vis, vis_transform, pan, pan_transform,
               pan_dtype, r_crs, dst_crs, weight,
               method="Brovey", src_nodata=0):
    """Pansharpen a lower-resolution visual band

    Parameters
    =========
    vis: ndarray, 3D with shape == (3, vh, vw)
        Visual band array with RGB bands
    vis_transform: Affine
        affine transform defining the georeferencing of the vis array
    pan: ndarray, 2D with shape == (ph, pw)
        Panchromatic band array
    pan_transform: Affine
        affine transform defining the georeferencing of the pan array
    method: string
        Algorithm for pansharpening; default Brovey

    Returns:
    ======
    pansharp: ndarray, 3D with shape == (3, ph, pw)
        pansharpened visual band
        affine transform is identical to `pan_transform`
    """
    rgb = _upsample(_create_apply_mask(vis), pan.shape, vis_transform, r_crs,
                    pan_transform, dst_crs)

    # Main Pansharpening Processing
    if method == "Brovey":
        pansharp, _ = Brovey(rgb, pan, weight, pan_dtype)
    # TODO: add other methods

    return pansharp 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def execute(self, eopatch):
        """ Execute method

        :param eopatch: input EOPatch
        :type eopatch: EOPatch
        :return: New EOPatch with vector data transformed into a raster feature
        :rtype: EOPatch
        """
        if eopatch.bbox is None:
            raise ValueError('EOPatch has to have a bounding box')

        height, width = self._get_raster_shape(eopatch)

        group_classes = self.overlap_value is not None
        rasterization_shapes = self._get_rasterization_shapes(eopatch, group_classes=group_classes)

        if not rasterization_shapes:
            eopatch[self.raster_feature] = np.full((height, width, 1), self.no_data_value, dtype=self.raster_dtype)
            return eopatch

        affine_transform = rasterio.transform.from_bounds(*eopatch.bbox, width=width, height=height)
        rasterize_args = dict(self.rasterio_params, transform=affine_transform, dtype=self.raster_dtype)

        raster = self._get_raster(eopatch, height, width)
        rasterize_func = rasterio.features.rasterize if self.overlap_value is None else self.rasterize_overlapped
        rasterize_func(rasterization_shapes, out=raster, **rasterize_args)

        eopatch[self.raster_feature] = raster[..., np.newaxis]
        return eopatch 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def _get_wms_request(self, bbox, size_x, size_y):
        """
        Returns WMS request.
        """
        bbox_3857 = bbox.transform(CRS.POP_WEB)

        return GeopediaWmsRequest(layer=self.layer,
                                  theme=self.theme,
                                  bbox=bbox_3857,
                                  width=size_x,
                                  height=size_y,
                                  image_format=self.image_format,
                                  custom_url_params={CustomUrlParam.TRANSPARENT: True}) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def accuracy(grid_in, guess_in, aoi_in, buffer_amount=0.01):
    """Measure accuracy against a specified grid 'truth' file.

    Parameters
    ----------
    grid_in : str, Path
        Path to vector truth file.
    guess_in : str, Path
        Path to guess output from guess2geom.
    aoi_in : str, Path
        Path to AOI feature.
    buffer_amount : float, optional (default 0.01.)
        Leeway in decimal degrees in calculating equivalence.
        0.01 DD equals approximately 1 mile at the equator.
    """

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

    grid = gpd.read_file(grid_in)
    grid_clipped = clip_line_poly(grid, aoi)
    grid_buff = grid_clipped.buffer(buffer_amount)

    guesses_reader = rasterio.open(guess_in)
    guesses = guesses_reader.read(1)

    grid_for_raster = [(row.geometry) for _, row in grid_clipped.iterrows()]
    grid_raster = rasterize(
        grid_for_raster,
        out_shape=guesses_reader.shape,
        fill=1,
        default_value=0,
        all_touched=True,
        transform=guesses_reader.transform,
    )
    grid_buff_raster = rasterize(
        grid_buff,
        out_shape=guesses_reader.shape,
        fill=1,
        default_value=0,
        all_touched=True,
        transform=guesses_reader.transform,
    )

    grid_raster = flip_arr_values(grid_raster)
    grid_buff_raster = flip_arr_values(grid_buff_raster)

    tp = true_positives(guesses, grid_buff_raster)
    fn = false_negatives(guesses, grid_raster)

    return tp, fn 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import transform [as 别名]
def _polygon_to_pix(polygon):
    """Transforms polygon coordinates to integer pixel coordinates. It makes
    the geometry easier to handle and reduces the number of points.

    Parameters
    ----------
    polygon: the shapely.geometry.Polygon instance to transform.

    Returns
    -------
    a shapely.geometry.Polygon class instance.
    """

    def project(x, y):
        return np.rint(x).astype(np.int64), np.rint(y).astype(np.int64)

    poly_pix = shapely.ops.transform(project, polygon)

    # simple trick to correct invalid polys:
    tmp = poly_pix.buffer(0)

    # sometimes the glacier gets cut out in parts
    if tmp.type == 'MultiPolygon':
        # If only small arms are cut out, remove them
        area = np.array([_tmp.area for _tmp in tmp])
        _tokeep = np.argmax(area).item()
        tmp = tmp[_tokeep]

        # check that the other parts really are small,
        # otherwise replace tmp with something better
        area = area / area[_tokeep]
        for _a in area:
            if _a != 1 and _a > 0.05:
                # these are extremely thin glaciers
                # eg. RGI40-11.01381 RGI40-11.01697 params.d1 = 5. and d2 = 8.
                # make them bigger until its ok
                for b in np.arange(0., 1., 0.01):
                    tmp = shapely.ops.transform(project, polygon.buffer(b))
                    tmp = tmp.buffer(0)
                    if tmp.type == 'MultiPolygon':
                        continue
                    if tmp.is_valid:
                        break
                if b == 0.99:
                    raise InvalidGeometryError('This glacier geometry is not '
                                               'valid.')

    if not tmp.is_valid:
        raise InvalidGeometryError('This glacier geometry is not valid.')

    return tmp