Python CRS.to_dict方法代码示例

# 需要导入模块: from rasterio.crs import CRS [as 别名]
# 或者: from rasterio.crs.CRS import to_dict [as 别名]
 def test_get_area_def_from_raster(self):
     from pyresample import utils
     from rasterio.crs import CRS
     from affine import Affine
     x_size = 791
     y_size = 718
     transform = Affine(300.0379266750948, 0.0, 101985.0,
                        0.0, -300.041782729805, 2826915.0)
     crs = CRS(init='epsg:3857')
     source = tmptiff(x_size, y_size, transform, crs=crs)
     area_id = 'area_id'
     proj_id = 'proj_id'
     name = 'name'
     area_def = utils._rasterio.get_area_def_from_raster(
         source, area_id=area_id, name=name, proj_id=proj_id)
     self.assertEqual(area_def.area_id, area_id)
     self.assertEqual(area_def.proj_id, proj_id)
     self.assertEqual(area_def.name, name)
     self.assertEqual(area_def.width, x_size)
     self.assertEqual(area_def.height, y_size)
     self.assertDictEqual(crs.to_dict(), area_def.proj_dict)
     self.assertTupleEqual(area_def.area_extent, (transform.c, transform.f + transform.e * y_size,
                                                  transform.c + transform.a * x_size, transform.f))

# 需要导入模块: from rasterio.crs import CRS [as 别名]
# 或者: from rasterio.crs.CRS import to_dict [as 别名]
def zones(
    input,
    output,
    variable,
    attribute,
    like,
    netcdf3,
    zip):

    """
    Create zones in a NetCDF from features in a shapefile.  This is intended
    to be used as input to zonal statistics functions; it is not intended
    as a direct replacement for rasterizing geometries into NetCDF.

    Only handles < 65,535 features for now.

    If --attribute is provided, any features that do not have this will not be
    assigned to zones.

    A values lookup will be used to store values.  The zones are indices of
    the unique values encountered when extracting features.
    The original values are stored in an additional variable with the name of
    the zones variable plus '_values'.

    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 = list(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:]
        out_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:
        if attribute:
            if not attribute in shp.meta['schema']['properties']:
                raise click.BadParameter('{0} not found in dataset'.format(attribute),
                                         param='--attribute', param_hint='--attribute')

            att_dtype = shp.meta['schema']['properties'][attribute].split(':')[0]
            if not att_dtype in ('int', 'str'):
                raise click.BadParameter('integer or string attribute required'.format(attribute),
                                         param='--attribute', param_hint='--attribute')

        transform_required = CRS(shp.crs) != template_crs
        geometries = []
        values = set()
        values_lookup = {}

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

        index = 0
        for f in shp.filter(bbox=bbox.as_list()):
            value = f['properties'].get(attribute) if attribute else int(f['id'])
            if value is not None:
                geom = f['geometry']
                if transform_required:
                    geom = transform_geom(shp.crs, template_crs, geom)

                geometries.append((geom, index))

                if not value in values:
                    values.add(value)
                    values_lookup[index] = value
                    index += 1

            # Otherwise, these will not be rasterized

        num_geometries = len(geometries)
        # Save a slot at the end for nodata
        if num_geometries < 255:
            dtype = numpy.dtype('uint8')
        elif num_geometries < 65535:
            dtype = numpy.dtype('uint16')
        else:
            raise click.UsageError('Too many features to rasterize: {0}, Exceptioning...'.format(num_geometries))

        fill_value = get_fill_value(dtype)
#.........这里部分代码省略.........

# 需要导入模块: from rasterio.crs import CRS [as 别名]
# 或者: from rasterio.crs.CRS import to_dict [as 别名]
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