Python warp.transform_bounds函数代码示例

def test_transform_bounds_densify():
    # This transform is non-linear along the edges, so densification produces
    # a different result than otherwise
    src_crs = {'init': 'EPSG:4326'}
    dst_crs = {'init': 'EPSG:32610'}
    assert numpy.allclose(
        transform_bounds(
            src_crs,
            dst_crs,
            -120, 40, -80, 64,
            densify_pts=0
        ),
        (
            646695.227266598, 4432069.056898901,
            4201818.984205882, 7807592.187464975
        )
    )

    assert numpy.allclose(
        transform_bounds(
            src_crs,
            dst_crs,
            -120, 40, -80, 64,
            densify_pts=100
        ),
        (
            646695.2272665979, 4432069.056898901,
            4201818.984205882, 7807592.187464977
        )
    )

def test_transform_bounds_densify_out_of_bounds():
    with pytest.raises(ValueError):
        transform_bounds(
            {'init': 'EPSG:4326'},
            {'init': 'EPSG:32610'},
            -120, 40, -80, 64,
            densify_pts=-10
        )

def test_transform_bounds():
    """CRSError is raised."""
    left, bottom, right, top = (
        -11740727.544603072, 4852834.0517692715, -11584184.510675032,
        5009377.085697309)
    src_crs = 'EPSG:3857'
    dst_crs = {'proj': 'foobar'}
    with pytest.raises(CRSError):
        transform_bounds(src_crs, dst_crs, left, bottom, right, top)

def test_transform_bounds_densify_out_of_bounds():
    with pytest.raises(ValueError):
        transform_bounds(
            {"init": "epsg:4326"},
            {"init": "epsg:32610"},
            -120,
            40,
            -80,
            64,
            densify_pts=-10,
        )

def test_transform_bounds_densify():
    # This transform is non-linear along the edges, so densification produces
    # a different result than otherwise
    src_crs = {"init": "epsg:4326"}
    dst_crs = {"init": "epsg:2163"}
    assert np.allclose(
        transform_bounds(src_crs, dst_crs, -120, 40, -80, 64, densify_pts=0),
        (-1684649.41338, -350356.81377, 1684649.41338, 2234551.18559),
    )

    assert np.allclose(
        transform_bounds(src_crs, dst_crs, -120, 40, -80, 64, densify_pts=100),
        (-1684649.41338, -555777.79210, 1684649.41338, 2234551.18559),
    )

def generate_chunk_tasks(image_source, tile_dim):
    tasks = []
    zoom = image_source.zoom
    (min_col, max_col) = (image_source.tile_bounds[0], image_source.tile_bounds[2])
    (min_row, max_row) = (image_source.tile_bounds[1], image_source.tile_bounds[3])

    for tile_col in range(min_col, min(max_col + 1, 2**zoom)):
        for tile_row in range(min_row, min(max_row + 1, 2**zoom)):
            tile_bounds = mercantile.bounds(tile_col, tile_row, zoom)
            (wm_left, wm_bottom, wm_right, wm_top)  = warp.transform_bounds("EPSG:4326",
                                                                           "EPSG:3857",
                                                                            tile_bounds.west,
                                                                            tile_bounds.south,
                                                                            tile_bounds.east,
                                                                            tile_bounds.north)
            affine = transform.from_bounds(wm_left, wm_bottom, wm_right, wm_top, tile_dim, tile_dim)
            target_meta = { 
                "transform": affine[:6],
                "width": tile_dim,
                "height": tile_dim 
            }

            target = os.path.join(image_source.image_folder, "%d/%d/%d.tif" % (zoom, tile_col, tile_row))
            task = ChunkTask(source_uri=image_source.source_uri,
                             target_meta=target_meta,
                             target=target)

            tasks.append(task)

    return tasks

def test_transform_bounds__esri_wkt():
    left, bottom, right, top = \
        (-78.95864996545055, 23.564991210854686,
         -76.57492370013823, 25.550873767433984)
    dst_projection_string = (
        'PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",'
        'GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",'
        'SPHEROID["GRS_1980",6378137.0,298.257222101]],'
        'PRIMEM["Greenwich",0.0],'
        'UNIT["Degree",0.0174532925199433]],'
        'PROJECTION["Albers"],'
        'PARAMETER["false_easting",0.0],'
        'PARAMETER["false_northing",0.0],'
        'PARAMETER["central_meridian",-96.0],'
        'PARAMETER["standard_parallel_1",29.5],'
        'PARAMETER["standard_parallel_2",45.5],'
        'PARAMETER["latitude_of_origin",23.0],'
        'UNIT["Meter",1.0],'
        'VERTCS["NAVD_1988",'
        'VDATUM["North_American_Vertical_Datum_1988"],'
        'PARAMETER["Vertical_Shift",0.0],'
        'PARAMETER["Direction",1.0],UNIT["Centimeter",0.01]]]')
    assert np.allclose(
        transform_bounds({"init": "epsg:4326"},
                         dst_projection_string,
                         left,
                         bottom,
                         right,
                         top),
        (
            1721263.7931814701,
            219684.49332178483,
            2002926.56696663,
            479360.16562217404),
    )

def test_transform_bounds_no_change():
    """Make sure that going from and to the same crs causes no change."""
    with rasterio.open('tests/data/RGB.byte.tif') as src:
        l, b, r, t = src.bounds
        assert np.allclose(
            transform_bounds(src.crs, src.crs, l, b, r, t),
            src.bounds
        )

def test_transform_bounds():
    with rasterio.open('tests/data/RGB.byte.tif') as src:
        l, b, r, t = src.bounds
        assert np.allclose(
            transform_bounds(src.crs, {'init': 'EPSG:4326'}, l, b, r, t),
            (
                -78.95864996545055, 23.564991210854686,
                -76.57492370013823, 25.550873767433984
            )
        )

def get_bounds(path):
    """ Retrun bounds in WGS84 system """

    with rasterio.drivers():
        src = rasterio.open(path)

        return transform_bounds(
                    src.crs,
                    {'init': 'EPSG:4326'},
                    *src.bounds)

    def clip(self):
        """ Clip images based on bounds provided
        Implementation is borrowed from
        https://github.com/brendan-ward/rasterio/blob/e3687ce0ccf8ad92844c16d913a6482d5142cf48/rasterio/rio/convert.py
        """

        self.output("Clipping", normal=True)

        # create new folder for clipped images
        path = check_create_folder(join(self.scene_path, 'clipped'))

        try:
            temp_bands = copy(self.bands)
            temp_bands.append('QA')
            for i, band in enumerate(temp_bands):
                band_name = self._get_full_filename(band)
                band_path = join(self.scene_path, band_name)

                self.output("Band %s" % band, normal=True, color='green', indent=1)
                with rasterio.open(band_path) as src:
                    bounds = transform_bounds(
                        {
                            'proj': 'longlat',
                            'ellps': 'WGS84',
                            'datum': 'WGS84',
                            'no_defs': True
                        },
                        src.crs,
                        *self.bounds
                    )

                    if disjoint_bounds(bounds, src.bounds):
                        bounds = adjust_bounding_box(src.bounds, bounds)

                    window = src.window(*bounds)

                    out_kwargs = src.meta.copy()
                    out_kwargs.update({
                        'driver': 'GTiff',
                        'height': window[0][1] - window[0][0],
                        'width': window[1][1] - window[1][0],
                        'transform': src.window_transform(window)
                    })

                    with rasterio.open(join(path, band_name), 'w', **out_kwargs) as out:
                        out.write(src.read(window=window))

            # Copy MTL to the clipped folder
            copyfile(join(self.scene_path, self.scene + '_MTL.txt'), join(path, self.scene + '_MTL.txt'))

            return path

        except IOError as e:
            exit(e.message, 1)

def test_transform_bounds():
    with rasterio.open("tests/data/RGB.byte.tif") as src:
        l, b, r, t = src.bounds
        assert np.allclose(
            transform_bounds(src.crs, {"init": "epsg:4326"}, l, b, r, t),
            (
                -78.95864996545055,
                23.564991210854686,
                -76.57492370013823,
                25.550873767433984,
            ),
        )

def _grid_datasets(datasets, bounds_override, grid_proj, grid_size):
    tiles = defaultdict(list)
    for dataset in datasets:
        dataset_proj = dataset.crs
        dataset_bounds = dataset.bounds
        bounds = bounds_override or BoundingBox(*transform_bounds(dataset_proj, grid_proj, *dataset_bounds))

        for y in range(int(bounds.bottom // grid_size[1]), int(bounds.top // grid_size[1]) + 1):
            for x in range(int(bounds.left // grid_size[0]), int(bounds.right // grid_size[0]) + 1):
                tile_index = (x, y)
                if _check_intersect(tile_index, grid_size, grid_proj, dataset_bounds, dataset_proj):
                    tiles[tile_index].append(dataset)

    return tiles

        def __call__(self):
            for i, path in enumerate(input):
                with rasterio.open(path) as src:
                    bounds = src.bounds
                    if dst_crs:
                        bbox = transform_bounds(src.crs,
                                                dst_crs, *bounds)
                    elif projection == 'mercator':
                        bbox = transform_bounds(src.crs,
                                                {'init': 'epsg:3857'}, *bounds)
                    elif projection == 'geographic':
                        bbox = transform_bounds(src.crs,
                                                {'init': 'epsg:4326'}, *bounds)
                    else:
                        bbox = bounds

                if precision >= 0:
                    bbox = [round(b, precision) for b in bbox]

                yield {
                    'type': 'Feature',
                    'bbox': bbox,
                    'geometry': {
                        'type': 'Polygon',
                        'coordinates': [[
                            [bbox[0], bbox[1]],
                            [bbox[2], bbox[1]],
                            [bbox[2], bbox[3]],
                            [bbox[0], bbox[3]],
                            [bbox[0], bbox[1]]]]},
                    'properties': {
                        'id': str(i),
                        'title': path,
                        'filename': os.path.basename(path)}}

                self._xs.extend(bbox[::2])
                self._ys.extend(bbox[1::2])

def clip_from_wgs(input_path, output_path, bounds):
    """
        Clips a raster based on WGS coordinates.
        Heavily borrowed from: https://github.com/mapbox/rasterio/blob/master/rasterio/rio/clip.py
        :param input_path The path to the raster to clip out of
        :param output_path The path to place the clipped raster
        :param bounds (left, bottom, right, top) bounds in WGS 84 projection
    """

    wgs_crs = CRS({'init': 'epsg:4326'})  # WGS 84 Projection

    with rasterio.open(input_path, 'r') as src:

        window = src.window(*transform_bounds(wgs_crs, src.crs, *bounds))
        height = window[0][1] - window[0][0]
        width = window[1][1] - window[1][0]
        t = 2048    # Threshold, for if the user selects an area greater than 2048 pixels in size
                    # If so, select the center of the area at size of 2048 by 2048
        if width > t:
            if width % 2 != 0:
                width -= 1
            trim = int((width - t) / 2)
            width = t
            window = (window[0], (window[1][0] + trim, window[1][0] + t + trim))

        if height > t:
            if height % 2 != 0:
                height -= 1
            trim = int((height - t) / 2)
            height = t
            window = ((window[0][0] + trim, window[0][0] + t + trim), window[1])

        out_kwargs = src.meta.copy()
        out_kwargs.update({
            'height': height,
            'width': width,
            'transform': src.window_transform(window),
            'dtype': 'int32',   # SyncroSim needs a signed int32
            'nodata': 0
            })

        with rasterio.open(output_path, 'w', **out_kwargs) as dst:
            dst.write(src.read(1, window=window).astype('int32'), 1)