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

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 all(round(x, 1) == x for x in flatten_coords(result["coordinates"]))

def mask(
    input,
    output,
    variable,
    like,
    netcdf3,
    all_touched,
    invert,
    zip):

    """
    Create a NetCDF mask from a shapefile.

    Values are equivalent to a numpy mask: 0 for unmasked areas, and 1 for masked areas.

    Template NetCDF dataset must have a valid projection defined or be inferred from dimensions (e.g., lat / long)
    """

    with Dataset(like) as template_ds:
        template_varname = data_variables(template_ds).keys()[0]
        template_variable = template_ds.variables[template_varname]
        template_crs = get_crs(template_ds, template_varname)

        if template_crs:
            template_crs = CRS.from_string(template_crs)
        elif is_geographic(template_ds, template_varname):
            template_crs = CRS({'init': 'EPSG:4326'})
        else:
            raise click.UsageError('template dataset must have a valid projection defined')

        spatial_dimensions = template_variable.dimensions[-2:]
        mask_shape = template_variable.shape[-2:]

        template_y_name, template_x_name = spatial_dimensions
        coords = SpatialCoordinateVariables.from_dataset(
            template_ds,
            x_name=template_x_name,
            y_name=template_y_name,
            projection=Proj(**template_crs.to_dict())
        )


    with fiona.open(input, 'r') as shp:
        transform_required = CRS(shp.crs) != template_crs

        # Project bbox for filtering
        bbox = coords.bbox
        if transform_required:
            bbox = bbox.project(Proj(**shp.crs), edge_points=21)

        geometries = []
        for f in shp.filter(bbox=bbox.as_list()):
            geom = f['geometry']
            if transform_required:
                geom = transform_geom(shp.crs, template_crs, geom)

            geometries.append(geom)

    click.echo('Converting {0} features to mask'.format(len(geometries)))

    if invert:
        fill_value = 0
        default_value = 1
    else:
        fill_value = 1
        default_value = 0

    with rasterio.Env():
        # Rasterize features to 0, leaving background as 1
        mask = rasterize(
            geometries,
            out_shape=mask_shape,
            transform=coords.affine,
            all_touched=all_touched,
            fill=fill_value,
            default_value=default_value,
            dtype=numpy.uint8
        )

    format = 'NETCDF3_CLASSIC' if netcdf3 else 'NETCDF4'
    dtype = 'int8' if netcdf3 else 'uint8'

    with Dataset(output, 'w', format=format) as out:
        coords.add_to_dataset(out, template_x_name, template_y_name)
        out_var = out.createVariable(variable, dtype, dimensions=spatial_dimensions, zlib=zip,
                                     fill_value=get_fill_value(dtype))
        out_var[:] = mask

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 all(round(x, 1) == x for x in flatten_coords(result['coordinates']))

def test_transform_geom_polygon_precision(polygon_3373):
    geom = polygon_3373
    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"]))

#.........这里部分代码省略.........

    precision: int (DEFAULT: -1)
        Decimal precision of coordinates. -1 for full float precision output

    Yields
    ------
    GeoJSON-like Feature dictionaries for shapes found in the given band
    """
    if bidx is not None and bidx > src.count:
        raise ValueError('bidx is out of range for raster')

    img = None
    msk = None

    # Adjust transforms.
    transform = src.transform
    if sampling > 1:
        # Determine the target shape (to decimate)
        shape = (int(math.ceil(src.height / sampling)),
                 int(math.ceil(src.width / sampling)))

        # Calculate independent sampling factors
        x_sampling = src.width / shape[1]
        y_sampling = src.height / shape[0]

        # Decimation of the raster produces a georeferencing
        # shift that we correct with a translation.
        transform *= Affine.translation(
            src.width % x_sampling, src.height % y_sampling)

        # And follow by scaling.
        transform *= Affine.scale(x_sampling, y_sampling)

    # Most of the time, we'll use the valid data mask.
    # We skip reading it if we're extracting every possible
    # feature (even invalid data features) from a band.
    if not band or (band and not as_mask and not with_nodata):
        if sampling == 1:
            msk = src.read_masks(bidx)
        else:
            msk_shape = shape
            if bidx is None:
                msk = np.zeros(
                    (src.count,) + msk_shape, 'uint8')
            else:
                msk = np.zeros(msk_shape, 'uint8')
            msk = src.read_masks(bidx, msk)

        if bidx is None:
            msk = np.logical_or.reduce(msk).astype('uint8')

        # Possibly overridden below.
        img = msk

    # Read the band data unless the --mask option is given.
    if band:
        if sampling == 1:
            img = src.read(bidx, masked=False)
        else:
            img = np.zeros(
                shape,
                dtype=src.dtypes[src.indexes.index(bidx)])
            img = src.read(bidx, img, masked=False)

    # If as_mask option was given, convert the image
    # to a binary image. This reduces the number of shape
    # categories to 2 and likely reduces the number of
    # shapes.
    if as_mask:
        tmp = np.ones_like(img, 'uint8') * 255
        tmp[img == 0] = 0
        img = tmp
        if not with_nodata:
            msk = tmp

    # Prepare keyword arguments for shapes().
    kwargs = {'transform': transform}
    if not with_nodata:
        kwargs['mask'] = msk

    src_basename = os.path.basename(src.name)

    # Yield GeoJSON features.
    for i, (g, val) in enumerate(
            rasterio.features.shapes(img, **kwargs)):
        if geographic:
            g = warp.transform_geom(
                src.crs, 'EPSG:4326', g,
                antimeridian_cutting=True, precision=precision)
        xs, ys = zip(*coords(g))
        yield {
            'type': 'Feature',
            'id': "{0}:{1}".format(src_basename, i),
            'properties': {
                'val': val,
                'filename': src_basename
            },
            'bbox': [min(xs), min(ys), max(xs), max(ys)],
            'geometry': g
        }