Python warp.transform_geom函数代码示例

def test_transform_geom_gdal22():
    """Enabling `antimeridian_cutting` has no effect on GDAL 2.2.0 or newer
    where antimeridian cutting is always enabled.  This could produce
    unexpected geometries, so an exception is raised.
    """
    geom = {"type": "Point", "coordinates": [0, 0]}
    with pytest.raises(GDALVersionError):
        transform_geom("EPSG:4326", "EPSG:3857", geom, antimeridian_cutting=False)

def test_transform_geom():
    geom = {
        'type': 'Polygon',
        'coordinates': (
            ((798842.3090855901, 6569056.500655151),
                (756688.2826828464, 6412397.888771972),
                (755571.0617232556, 6408461.009397383),
                (677605.2284582685, 6425600.39266733),
                (677605.2284582683, 6425600.392667332),
                (670873.3791649605, 6427248.603432341),
                (664882.1106069803, 6407585.48425362),
                (663675.8662823177, 6403676.990080649),
                (485120.71963574126, 6449787.167760638),
                (485065.55660851026, 6449802.826920689),
                (485957.03982722526, 6452708.625101285),
                (487541.24541826674, 6457883.292107048),
                (531008.5797472061, 6605816.560367976),
                (530943.7197027118, 6605834.9333479265),
                (531888.5010308184, 6608940.750411527),
                (533299.5981959199, 6613962.642851984),
                (533403.6388841148, 6613933.172096095),
                (576345.6064638699, 6761983.708069147),
                (577649.6721159086, 6766698.137844516),
                (578600.3589008929, 6770143.99782289),
                (578679.4732294685, 6770121.638265098),
                (655836.640492081, 6749376.357102599),
                (659913.0791150068, 6764770.1314677475),
                (661105.8478791204, 6769515.168134831),
                (661929.4670843681, 6772800.8565198565),
                (661929.4670843673, 6772800.856519875),
                (661975.1582566603, 6772983.354777632),
                (662054.7979028501, 6772962.86384242),
                (841909.6014891531, 6731793.200435557),
                (840726.455490463, 6727039.8672589315),
                (798842.3090855901, 6569056.500655151)),
            )
    }

    result = transform_geom('EPSG:3373', 'EPSG:4326', geom)
    assert result['type'] == 'Polygon'
    assert len(result['coordinates']) == 1

    result = transform_geom(
        'EPSG:3373', 'EPSG:4326', geom, antimeridian_cutting=True)
    assert result['type'] == 'MultiPolygon'
    assert len(result['coordinates']) == 2

    result = transform_geom(
        'EPSG:3373', 
        'EPSG:4326', 
        geom, 
        antimeridian_cutting=True, 
        antimeridian_offset=0)
    assert result['type'] == 'MultiPolygon'
    assert len(result['coordinates']) == 2

    result = transform_geom('EPSG:3373', 'EPSG:4326',  geom,  precision=1)
    assert int(result['coordinates'][0][0][0] * 10) == -1778

def test_transform_geom_gdal22():
    """Enabling `antimeridian_cutting` has no effect on GDAL 2.2.0 or newer
    where antimeridian cutting is always enabled.  This could produce
    unexpected geometries, so an exception is raised.
    """
    geom = {
        'type': 'Point',
        'coordinates': [0, 0]
    }
    with pytest.raises(GDALBehaviorChangeException):
        transform_geom(
            'EPSG:4326', 'EPSG:3857', geom, antimeridian_cutting=False)

def test_transform_geom_polygon_offset(polygon_3373):
    geom = polygon_3373
    result = transform_geom(
        "EPSG:3373", "EPSG:4326", geom, antimeridian_cutting=True, antimeridian_offset=0
    )
    assert result["type"] == "MultiPolygon"
    assert len(result["coordinates"]) == 2

def test_transform_geom_linearring_precision(polygon_3373):
    ring = polygon_3373["coordinates"][0]
    geom = {"type": "LinearRing", "coordinates": ring}
    result = transform_geom(
        "EPSG:3373", "EPSG:4326", geom, precision=1, antimeridian_cutting=True
    )
    assert all(round(x, 1) == x for x in flatten_coords(result["coordinates"]))

def test_transform_geom_linestring_precision_z(polygon_3373):
    ring = polygon_3373["coordinates"][0]
    x, y = zip(*ring)
    ring = list(zip(x, y, [0.0 for i in range(len(x))]))
    geom = {"type": "LineString", "coordinates": ring}
    result = transform_geom("EPSG:3373", "EPSG:3373", geom, precision=1)
    assert int(result["coordinates"][0][0] * 10) == 7988423
    assert int(result["coordinates"][0][2] * 10) == 0

def test_issue_1446():
    """Confirm resolution of #1446"""
    g = transform_geom(
        CRS.from_epsg(4326),
        CRS.from_epsg(32610),
        {"type": "Point", "coordinates": (-122.51403808499907, 38.06106733107932)},
    )
    assert round(g["coordinates"][0], 1) == 542630.9
    assert round(g["coordinates"][1], 1) == 4212702.1

def test_transform_geom_polygon_offset(polygon_3373):
    geom = polygon_3373
    result = transform_geom(
        'EPSG:3373',
        'EPSG:4326',
        geom,
        antimeridian_cutting=True,
        antimeridian_offset=0)
    assert result['type'] == 'MultiPolygon'
    assert len(result['coordinates']) == 2

    def _add_zones(self, source, name, uid, species, path, zone_field):
        with fiona.open(path, 'r') as shp:
            for feature in shp:
                try:
                    zone_id = feature['properties'][zone_field]
                except KeyError:
                    zone_id = feature['properties'][zone_field.lower()]

                if zone_id is None:
                    continue

                zone_id = int(zone_id)

                if zone_id == 0:
                    continue

                if hasattr(name, '__call__'):
                    object_id = feature['properties'].get('OBJECTID')
                    zone_name = name(zone_id, object_id)
                else:
                    zone_name = name.format(zone_id=zone_id, species=SPECIES_NAMES.get(species))

                geometry = transform_geom(shp.crs, {'init': 'EPSG:4326'}, feature['geometry'])

                if feature['geometry']['type'] == 'MultiPolygon':
                    geometry['coordinates'] = itertools.chain(*geometry['coordinates'])

                polygon = Polygon(*[LinearRing(x) for x in geometry['coordinates']])
                uid_suffix = 0

                while True:
                    zone_uid = uid.format(zone_id=zone_id)

                    if uid_suffix > 0:
                        zone_uid += '_{}'.format(uid_suffix)

                    try:
                        with transaction.atomic():
                            SeedZone.objects.create(
                                source=source, name=zone_name, species=species, zone_id=zone_id, zone_uid=zone_uid,
                                polygon=polygon
                            )
                        break
                    except IntegrityError:
                        if uid_suffix > 10:
                            raise

                        uid_suffix += 1

    def handle(self, name, file, *args, **options):
        name = name[0]
        file = file[0]

        if Region.objects.filter(name__iexact=name).exists():
            message = (
                'WARNING: This will replace an existing region with the same name: {}. Do you want to continue? [y/n]'
            ).format(name)
            if input(message).lower() not in {'y', 'yes'}:
                return

        temp_dir = None

        try:
            if file.endswith('.zip'):
                temp_dir = mkdtemp()

                with ZipFile(file) as zf:
                    zf.extractall(temp_dir)

                    try:
                        file = glob.glob(os.path.join(temp_dir, '*.shp'))[0]
                    except IndexError:
                        raise ValueError('No shapefile in zip archive')

            polygons = []

            with fiona.open(file, 'r') as shp:
                for feature in shp:
                    geometry = transform_geom(shp.crs, {'init': 'EPSG:4326'}, feature['geometry'])
                    polygons.append(Polygon(*[LinearRing(x) for x in geometry['coordinates']]))

            with transaction.atomic():
                Region.objects.filter(name__iexact=name).delete()
                Region.objects.create(name=name, polygons=MultiPolygon(polygons))

        finally:
            if temp_dir is not None:
                try:
                    shutil.rmtree(temp_dir)
                except OSError:
                    pass

def test_issue_1446_b():
    """Confirm that lines aren't thrown as reported in #1446"""
    src_crs = CRS({"init": "epsg:4326"})
    dst_crs = CRS(
        {
            "proj": "sinu",
            "lon_0": 350.85607029556,
            "x_0": 0,
            "y_0": 0,
            "a": 3396190,
            "b": 3396190,
            "units": "m",
            "no_defs": True,
        }
    )
    collection = json.load(open("tests/data/issue1446.geojson"))
    geoms = {f["properties"]["fid"]: f["geometry"] for f in collection["features"]}
    transformed_geoms = {
        k: transform_geom(src_crs, dst_crs, g) for k, g in geoms.items()
    }
    # Before the fix, this geometry was thrown eastward of 0.0. It should be between -350 and -250.
    assert all([-350 < x < -150 for x, y in transformed_geoms[183519]["coordinates"]])

def test_transform_geom_linestring_precision(polygon_3373):
    ring = polygon_3373['coordinates'][0]
    geom = {'type': 'LineString', 'coordinates': ring}
    result = transform_geom('EPSG:3373', 'EPSG:4326', geom, precision=1, antimeridian_cutting=True)
    assert all(round(x, 1) == x for x in flatten_coords(result['coordinates']))

def test_transform_geom_dst_crs_none():
    with pytest.raises(CRSError):
        transform_geom(WGS84_crs, None, None)

def test_transform_geom_multipolygon(polygon_3373):
    geom = {
        'type': 'MultiPolygon', 'coordinates': [polygon_3373['coordinates']]}
    result = transform_geom('EPSG:3373', 'EPSG:4326', geom, precision=1)
    assert int(result['coordinates'][0][0][0][0] * 10) == -1778

def test_transform_geom_polygon_precision(polygon_3373):
    geom = polygon_3373
    result = transform_geom('EPSG:3373', 'EPSG:4326', geom, precision=1)
    assert int(result['coordinates'][0][0][0] * 10) == -1778