Python fiona.open() Examples

def ConvertRaster2LatLong(InputRasterFile,OutputRasterFile):

    """
    Convert a raster to lat long WGS1984 EPSG:4326 coordinates for global plotting

    MDH

    """

    # import modules
    import rasterio
    from rasterio.warp import reproject, calculate_default_transform as cdt, Resampling

    # read the source raster
    with rasterio.open(InputRasterFile) as src:
        #get input coordinate system
        Input_CRS = src.crs
        # define the output coordinate system
        Output_CRS = {'init': "epsg:4326"}
        # set up the transform
        Affine, Width, Height = cdt(Input_CRS,Output_CRS,src.width,src.height,*src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({
            'crs': Output_CRS,
            'transform': Affine,
            'affine': Affine,
            'width': Width,
            'height': Height
        })

        with rasterio.open(OutputRasterFile, 'w', **kwargs) as dst:
            for i in range(1, src.count+1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.affine,
                    src_crs=src.crs,
                    dst_transform=Affine,
                    dst_crs=Output_CRS,
                    resampling=Resampling.bilinear) 

def wkt_to_shp(wkt_list, shp_file):
    '''Take output of build_graph_wkt() and render the list of linestrings
    into a shapefile
    # https://gis.stackexchange.com/questions/52705/how-to-write-shapely-geometries-to-shapefiles
    '''

    # Define a linestring feature geometry with one attribute
    schema = {
        'geometry': 'LineString',
        'properties': {'id': 'int'},
    }

    # Write a new shapefile
    with fiona.open(shp_file, 'w', 'ESRI Shapefile', schema) as c:
        for i, line in enumerate(wkt_list):
            shape = shapely.wkt.loads(line)
            c.write({
                    'geometry': mapping(shape),
                    'properties': {'id': i},
                    })

    return


############################################################################### 

def polygonize_raster(infile, outfile, mask_value=1, driver='ESRI Shapefile'):

    with rasterio.open(infile) as src:

        image = src.read(1)

        if mask_value is not None:
            mask = image == mask_value
        else:
            mask = None

        results = (
            {'properties': {'raster_val': v}, 'geometry': s}
            for i, (s, v)
            in enumerate(
                shapes(image, mask=mask, transform=src.transform)))

        with fiona.open(
            outfile, 'w',
            driver=driver,
            crs=src.crs,
            schema={'properties': [('raster_val', 'int')],
                    'geometry': 'Polygon'}) as dst:

            dst.writerecords(results) 

def test_convert_single_gtiff(cleantopo_br_tif, mp_tmpdir):
    """Automatic geodetic tile pyramid creation of raster files."""
    single_gtiff = os.path.join(mp_tmpdir, "single_out.tif")
    run_cli([
        "convert",
        cleantopo_br_tif,
        single_gtiff,
        "--output-pyramid", "geodetic",
        "-z", "3"
    ])
    with rasterio.open(single_gtiff, "r") as src:
        assert src.meta["driver"] == "GTiff"
        assert src.meta["dtype"] == "uint16"
        data = src.read(masked=True)
        assert data.mask.any()
        assert not src.overviews(1) 

def test_convert_single_gtiff_overviews(cleantopo_br_tif, mp_tmpdir):
    """Automatic geodetic tile pyramid creation of raster files."""
    single_gtiff = os.path.join(mp_tmpdir, "single_out.tif")
    run_cli([
        "convert",
        cleantopo_br_tif,
        single_gtiff,
        "--output-pyramid", "geodetic",
        "-z", "3",
        "--overviews",
        "--overviews-resampling-method", "bilinear"
    ])
    with rasterio.open(single_gtiff, "r") as src:
        assert src.meta["driver"] == "GTiff"
        assert src.meta["dtype"] == "uint16"
        data = src.read(masked=True)
        assert data.mask.any()
        assert src.overviews(1) 

def read_terrace_centrelines(DataDirectory, shapefile_name):
    """
    This function reads in a shapefile of terrace centrelines
    using shapely and fiona

    Args:
        DataDirectory (str): the data directory
        shapefile_name (str): the name of the shapefile

    Returns: shapely polygons with terraces

    Author: FJC
    """
    Lines = {}
    with fiona.open(DataDirectory+shapefile_name, 'r') as input:
        for f in input:
            this_line = LineString(shape(f['geometry']))
            this_id = f['properties']['id']
            Lines[this_id] = this_line
    return Lines 

def read_terrace_shapefile(DataDirectory, shapefile_name):
    """
    This function reads in a shapefile of digitised terraces
    using shapely and fiona

    Args:
        DataDirectory (str): the data directory
        shapefile_name (str): the name of the shapefile

    Returns: shapely polygons with terraces

    Author: FJC
    """
    Polygons = {}
    with fiona.open(DataDirectory+shapefile_name, 'r') as input:
        for f in input:
            this_shape = Polygon(shape(f['geometry']))
            this_id = f['properties']['id']
            Polygons[this_id] = this_shape

    return Polygons 

def test_convert_dtype(cleantopo_br_tif, mp_tmpdir):
    """Automatic tile pyramid creation using dtype scale."""
    run_cli([
        "convert",
        cleantopo_br_tif,
        mp_tmpdir,
        "--output-pyramid", "mercator",
        "--output-dtype", "uint8"
    ])
    for zoom, row, col in [(4, 15, 15), (3, 7, 7)]:
        out_file = os.path.join(*[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
        with rasterio.open(out_file, "r") as src:
            assert src.meta["driver"] == "GTiff"
            assert src.meta["dtype"] == "uint8"
            data = src.read(masked=True)
            assert data.mask.any() 

def test_convert_scale_ratio(cleantopo_br_tif, mp_tmpdir):
    """Automatic tile pyramid creation cropping data."""
    run_cli([
        "convert",
        cleantopo_br_tif,
        mp_tmpdir,
        "--output-pyramid", "mercator",
        "--output-dtype", "uint8",
        "--scale-ratio", "0.003"
    ])
    for zoom, row, col in [(4, 15, 15), (3, 7, 7)]:
        out_file = os.path.join(*[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
        with rasterio.open(out_file, "r") as src:
            assert src.meta["driver"] == "GTiff"
            assert src.meta["dtype"] == "uint8"
            data = src.read(masked=True)
            assert data.mask.any()
            assert not data.mask.all() 

def test_convert_scale_offset(cleantopo_br_tif, mp_tmpdir):
    """Automatic tile pyramid creation cropping data."""
    run_cli([
        "convert",
        cleantopo_br_tif,
        mp_tmpdir,
        "--output-pyramid", "mercator",
        "--output-dtype", "uint8",
        "--scale-offset", "1"
    ])
    for zoom, row, col in [(4, 15, 15), (3, 7, 7)]:
        out_file = os.path.join(*[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
        with rasterio.open(out_file, "r") as src:
            assert src.meta["driver"] == "GTiff"
            assert src.meta["dtype"] == "uint8"
            data = src.read(masked=True)
            assert data.mask.any()
            assert not data.mask.all() 

def test_convert_mapchete(cleantopo_br, mp_tmpdir):
    # prepare data
    with mapchete.open(cleantopo_br.path) as mp:
        mp.batch_process(zoom=[1, 4])
    run_cli([
        "convert",
        cleantopo_br.path,
        mp_tmpdir,
        "--output-pyramid", "geodetic",
        "--output-metatiling", "1",
        "-d",
    ])
    for zoom, row, col in [(4, 15, 31), (3, 7, 15), (2, 3, 7), (1, 1, 3)]:
        out_file = os.path.join(*[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
        with rasterio.open(out_file, "r") as src:
            assert src.meta["driver"] == "GTiff"
            assert src.meta["dtype"] == "uint16"
            data = src.read(masked=True)
            assert data.mask.any() 

def test_convert_tiledir(cleantopo_br, mp_tmpdir):
    # prepare data
    with mapchete.open(cleantopo_br.path) as mp:
        mp.batch_process(zoom=[1, 4])
    run_cli([
        "convert",
        os.path.join(
            cleantopo_br.dict["config_dir"], cleantopo_br.dict["output"]["path"]
        ),
        mp_tmpdir,
        "--output-pyramid", "geodetic",
        "--output-metatiling", "1",
        "--zoom", "1,4",
        "-d",
    ])
    for zoom, row, col in [(4, 15, 31), (3, 7, 15), (2, 3, 7), (1, 1, 3)]:
        out_file = os.path.join(*[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
        with rasterio.open(out_file, "r") as src:
            assert src.meta["driver"] == "GTiff"
            assert src.meta["dtype"] == "uint16"
            data = src.read(masked=True)
            assert data.mask.any() 

def test_convert_single_gtiff_cog(cleantopo_br_tif, mp_tmpdir):
    """Automatic geodetic tile pyramid creation of raster files."""
    single_gtiff = os.path.join(mp_tmpdir, "single_out_cog.tif")
    run_cli([
        "convert",
        cleantopo_br_tif,
        single_gtiff,
        "--output-pyramid", "geodetic",
        "-z", "3",
        "--cog"
    ])
    with rasterio.open(single_gtiff, "r") as src:
        assert src.meta["driver"] == "GTiff"
        assert src.meta["dtype"] == "uint16"
        data = src.read(masked=True)
        assert data.mask.any()
    assert cog_validate(single_gtiff, strict=True) 

def test_convert_png(cleantopo_br_tif, mp_tmpdir):
    """Automatic PNG tile pyramid creation of raster files."""
    run_cli([
        "convert",
        cleantopo_br_tif,
        mp_tmpdir,
        "--output-pyramid", "mercator",
        "--output-format", "PNG"
    ])
    for zoom, row, col in [(4, 15, 15), (3, 7, 7)]:
        out_file = os.path.join(
            *[mp_tmpdir, str(zoom), str(row), str(col) + ".png"])
        with rasterio.open(out_file, "r") as src:
            assert src.meta["driver"] == "PNG"
            assert src.meta["dtype"] == "uint8"
            data = src.read(masked=True)
            assert data.mask.any() 

def upload_s3(bucket, json_file, metadata):
    """
    파일을 gz하여 s3로 업로드
    :param json_file: 업로드할 json 파일명
    :return:
    """
    gz_name = f"{json_file}.gz"
    obj_key = f"json/{path.basename(gz_name)}"
    print("업로드", gz_name, obj_key)

    with open(json_file, 'rb') as f:
        gz = gzip.compress(f.read())
        s3.put_object(
            Body=gz,
            Bucket=bucket,
            ContentEncoding='gzip',
            ContentLanguage='string',
            ContentType='application/json',
            Key=obj_key,
            # todo : 메타데이터 추가 - 2018-07-28
            Metadata=metadata,
        ) 

def read(self, output_tile, **kwargs):
        """
        Read existing process output.

        Parameters
        ----------
        output_tile : ``BufferedTile``
            must be member of output ``TilePyramid``

        Returns
        -------
        process output : list
        """
        try:
            with fiona.open(self.get_path(output_tile), "r") as src:
                return list(src)
        except DriverError as e:
            for i in ("does not exist in the file system", "No such file or directory"):
                if i in str(e):
                    return self.empty(output_tile)
            else:
                raise 

def bbox(self, out_crs=None):
        """
        Return data bounding box.

        Parameters
        ----------
        out_crs : ``rasterio.crs.CRS``
            rasterio CRS object (default: CRS of process pyramid)

        Returns
        -------
        bounding box : geometry
            Shapely geometry object
        """
        out_crs = self.pyramid.crs if out_crs is None else out_crs
        with fiona.open(self.path) as inp:
            inp_crs = CRS(inp.crs)
            bbox = box(*inp.bounds)
        # TODO find a way to get a good segmentize value in bbox source CRS
        return reproject_geometry(bbox, src_crs=inp_crs, dst_crs=out_crs) 

def test_shp_to_geojson_records_geom_type_is_multilinestring(
    create_input_file, create_output_centerline_file
):
    EXPECTED_TYPE = "MultiLineString"

    input_polygon_shp = create_input_file("polygons", "shp")
    output_centerline_geojson = create_output_centerline_file("geojson")

    runner = CliRunner()
    runner.invoke(
        create_centerlines, [input_polygon_shp, output_centerline_geojson]
    )

    with fiona.open(output_centerline_geojson) as dst:
        for record in dst:
            assert record.get("geometry").get("type") == EXPECTED_TYPE 

def test_shp_to_shp_records_geom_type_is_multilinestring(
    create_input_file, create_output_centerline_file
):
    EXPECTED_TYPE = "MultiLineString"

    input_polygon_shp = create_input_file("polygons", "shp")
    output_centerline_shp = create_output_centerline_file("shp")

    runner = CliRunner()
    runner.invoke(
        create_centerlines, [input_polygon_shp, output_centerline_shp]
    )

    with fiona.open(output_centerline_shp) as dst:
        for record in dst:
            assert record.get("geometry").get("type") == EXPECTED_TYPE 

def overlay_bitmap(bitmap, raster_dataset, out_path, color='blue'):
    """Overlay the given satellite image with a bitmap."""

    # RGB values for possible color options.
    colors = {
        "red": (255, 0, 0),
        "green": (0, 255, 0),
        "blue": (0, 0, 255)
    }

    red, green, blue = raster_dataset.read()
    red[bitmap == 1] = colors[color][0]
    green[bitmap == 1] = colors[color][1]
    blue[bitmap == 1] = colors[color][2]

    profile = raster_dataset.profile
    with rasterio.open(out_path, 'w', **profile) as dst:
        dst.write(red, 1)
        dst.write(green, 2)
        dst.write(blue, 3)

    return rasterio.open(out_path) 

def test_read_raster_window_resampling(cleantopo_br_tif):
    """Assert various resampling options work."""
    tp = BufferedTilePyramid("geodetic")
    with rasterio.open(cleantopo_br_tif, "r") as src:
        tiles = tp.tiles_from_bounds(src.bounds, 4)
    for tile in tiles:
        outputs = [
            read_raster_window(
                cleantopo_br_tif, tile, resampling=resampling)
            for resampling in [
                "nearest", "bilinear", "cubic", "cubic_spline", "lanczos",
                "average", "mode"
            ]
        ]
        # resampling test:
        assert any([
            not np.array_equal(w, v) for v, w in zip(outputs[:-1], outputs[1:])
        ]) 

def shp2meta(s3_obj):
    tmp_file = path.join(tmp, s3_obj.object_key)
    file_name = s3_obj.object_key.split('.')[0]
    # 파일 다운로드
    with open(tmp_file, 'wb') as data:
        s3_obj.object.download_fileobj(data)

    # 압축 해제
    extract_dir = path.join(tmp, file_name)
    make_dir(extract_dir)
    with tarfile.open(tmp_file) as tar:
        tar.extractall(extract_dir)

    # json 파일로 변경
    shp_file = [path.join(extract_dir, name) for name in os.listdir(extract_dir) if name.split('.')[-1] == "shp"][0]
    metadata = get_meta(shp_file)
    s3_meta = s3_obj.object.get()['Metadata']
    filter_keys = ['upload_by', 'uploader', 'description']
    for k in filter_keys:
        for key in s3_meta:
            if k in key:
                metadata[k] = s3_meta[key]
                continue
    record = MetaData(hash_key=s3_obj.object_key, **metadata)
    record.save() 

def test_convert_bidx(cleantopo_br_tif, mp_tmpdir):
    """Automatic geodetic tile pyramid creation of raster files."""
    single_gtiff = os.path.join(mp_tmpdir, "single_out_bidx.tif")
    run_cli([
        "convert",
        cleantopo_br_tif,
        single_gtiff,
        "--output-pyramid", "geodetic",
        "-z", "3",
        "--bidx", "1"
    ])
    with rasterio.open(single_gtiff, "r") as src:
        assert src.meta["driver"] == "GTiff"
        assert src.meta["dtype"] == "uint16"
        data = src.read(masked=True)
        assert data.mask.any()
        assert not src.overviews(1) 

def parse_gpx(self, filename):
        tracks_layer = fiona.open(filename, layer='tracks')
        feature      = tracks_layer[0]
        self.points  = feature['geometry']['coordinates']
        self.track   = self.check_track(MultiLineString(self.points))
        self.name    = feature['properties']['name'] 

def open(self, tile, **kwargs):
        """
        Return InputTile object.

        Parameters
        ----------
        tile : ``Tile``

        Returns
        -------
        input tile : ``InputTile``
            tile view of input data
        """
        return InputTile(tile, self, **kwargs) 

def open(self, tile, process):
        """
        Open process output as input for other process.

        Parameters
        ----------
        tile : ``Tile``
        process : ``MapcheteProcess``
        """
        return InputTile(tile, process) 

def test_convert_geodetic(cleantopo_br_tif, mp_tmpdir):
    """Automatic geodetic tile pyramid creation of raster files."""
    run_cli(["convert", cleantopo_br_tif, mp_tmpdir, "--output-pyramid", "geodetic"])
    for zoom, row, col in [(4, 15, 31), (3, 7, 15), (2, 3, 7), (1, 1, 3)]:
        out_file = os.path.join(
            *[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
        with rasterio.open(out_file, "r") as src:
            assert src.meta["driver"] == "GTiff"
            assert src.meta["dtype"] == "uint16"
            data = src.read(masked=True)
            assert data.mask.any() 

def driver_from_file(input_file):
    """
    Guess driver from file extension.

    Returns
    -------
    driver : string
        driver name
    """
    file_ext = os.path.splitext(input_file)[1].split(".")[1]
    if file_ext == "mapchete":
        return "Mapchete"
    try:
        with rasterio.open(input_file):
            return "raster_file"
    except:
        try:
            with fiona.open(input_file):
                return "vector_file"
        except:
            if path_exists(input_file):
                raise MapcheteDriverError(
                    "%s has an unknown file extension or could not be opened by neither "
                    "rasterio nor fiona." % input_file
                )
            else:
                raise FileNotFoundError("%s does not exist" % input_file) 

def vertexcount(input):
    with open(input) as aoi:
        aoi_resp = json.load(aoi)
        if list_depth(aoi_resp['features'][0]['geometry']['coordinates'])==0:
            print('Number of vertices in simplified geometry: '+str(len(aoi_resp['features'][0]['geometry']['coordinates'][0])))
        elif list_depth(aoi_resp['features'][0]['geometry']['coordinates'])==1:
            print('Number of vertices in simplfied geometry: '+str(len(aoi_resp['features'][0]['geometry']['coordinates'][0][0])))
        else:
            print('Please check GeoJSON: Could not parse coordinates') 

def test_raster_window_memoryfile():
    """Use context manager for rasterio MemoryFile."""
    tp = BufferedTilePyramid("geodetic")
    tile = tp.tile(5, 5, 5)
    data = ma.masked_array(np.ones((2, ) + tile.shape))
    for out_profile in [
        dict(
            driver="GTiff", count=2, dtype="uint8", compress="lzw", nodata=0,
            height=tile.height, width=tile.width, affine=tile.affine),
        dict(
            driver="GTiff", count=2, dtype="uint8", compress="deflate",
            nodata=0, height=tile.height, width=tile.width,
            affine=tile.affine),
        dict(
            driver="PNG", count=2, dtype="uint8", nodata=0, height=tile.height,
            width=tile.width, compress=None, affine=tile.affine),
    ]:
        with RasterWindowMemoryFile(
            in_tile=tile, in_data=data, out_profile=out_profile
        ) as memfile:
            with memfile.open() as src:
                assert src.read().any()
                assert src.meta["driver"] == out_profile["driver"]
                assert src.transform == tile.affine
                if out_profile["compress"]:
                    assert src.compression == Compression(
                        out_profile["compress"].upper())