[GIS] Take pixel locations from rasterio raster and write out shapefile

Question

I'm pulling a satellite image into python with rasterio, running some analysis to find a specific bounding box within that image, and then want to write those points out as a shapefile. I'm struggling with the Fiona documentation, much online seems out of date and I can't convert the official docs well to my example.

So basically I've got

import rasterio

dataset = rasterio.open(image_filepath)
print(dataset.crs)

# this outputs '+init=epsg:4326'

"""
run analysis to find points on raster
"""

# example of a single rectangle that we get from that analysis
geolocated_xs = (-113.21645410141903, -113.21645410141903, -113.2167506315605, -113.2167506315605) 
geolocated_ys = (26.827766731681987, 26.828098563030775, 26.827766731681987, 26.828098563030775)

"""
write out to shapefile
I'm not totally sure what this should look like...
"""

# I assume I need shapely somewhere here?
from shapely.geometry import mapping, Polygon

import fiona
from fiona.crs import from_epsg
crs = from_epsg(4326)

schema = ?

# so from the docs

with fiona.open(
     'foo.shp',
     'w',
     driver='ESRI Shapefile,
     crs=crs,
     schema=) as c:
        c.write( ??? )

But I'm not sure what to write out and how I need to convert those tuples to have them go into this Shapefile.

I seek suggestions for pointers to good examples that I can use to learn.

Answer 1

You're almost there. Updated docs are at https://fiona.readthedocs.io

Your schema is just a simple dict with 'geometry' and 'properties' (for the attribute fields) items.

You then need to populate a record with the coordinates as x,y tuples (no shapely required for something simple like this) and whatever attributes you like.

The coordinate tuples/lists for a polygon need to be nested:

• polygon
  • ring/s i.e holes (not sure how multipolygons are represented)
    • points

So a geometry would look something like:

[((x,y),(x,y),(x,y),(x,y))] # Note use of list for single element, can use a tuple, but that requires a trailing comma (((x,y),(x,y etc)),)

Just be careful with coordinate ordering as I've noted in the snippet below.

from collections import OrderedDict
import fiona
from fiona.crs import from_epsg

# example of a single rectangle that we get from that analysis
geolocated_xs = (-113.21645410141903, -113.21645410141903, -113.2167506315605, -113.2167506315605)
geolocated_ys = (26.827766731681987, 26.828098563030775, 26.827766731681987, 26.828098563030775)

# You could just zip the geolocated_xys together to form your polygon, 
# i.e. polygon = [zip(geolocated_xs, geolocated_ys)]

# but you need to be careful with ordering as the ordering of your
# coordinates will give a bow-tie self intersection
# (they are currently in lr->ur->ll->ul order)

# Fix the ordering otherwise you'll get a bow-tie self intersection
lr, ur, ll, ul = zip(geolocated_xs, geolocated_ys)
polygon = [(lr, ur, ul, ll)]

crs = from_epsg(4326)

# Define your schema as a polygon geom with a couple of fields
schema = {
    'geometry': 'Polygon',
    'properties': OrderedDict([
        ('somecol', 'str'),
        ('anothercol', 'int')
  ])
}

with fiona.open(
    'oo.shp',
    'w',
    driver='ESRI Shapefile',
    crs=crs,
    schema=schema) as c:

    record = {
        'geometry': {'coordinates': polygon, 'type': 'Polygon'},
        'id': '0',
        'properties': OrderedDict([('somecol', 'Some Value'),
                                   ('anothercol', 12345)
                                   ]),
        'type': 'Feature'}

    c.write(record)

And for multiple records:

# example of a multiple rectangles using your coordinate ordering
geolocated_xys = (
    (
        (-113.21645410141903, -113.21645410141903, -113.2167506315605, -113.2167506315605),
        (26.827766731681987, 26.828098563030775, 26.827766731681987, 26.828098563030775)
    ),
    (
        (-113.21615757127756, -113.21615757127756, -113.21645410141903, -113.21645410141903),
        (26.827766731681987, 26.828098563030775, 26.827766731681987, 26.828098563030775)
    )
)

etc...

with fiona.open(
    'c:/temp/oo.shp',
    'w',
    driver='ESRI Shapefile',
    crs=crs,
    schema=schema) as c:

    for id, polygon in enumerate(geolocated_xys):
        # Fix the ordering otherwise you'll get a bow-tie self intersection
        lr, ur, ll, ul = zip(*polygon)
        polygon = [(lr, ur, ul, ll)]

        record = {
            'geometry': {'coordinates': polygon, 'type': 'Polygon'},
            'id': id,
            'properties': OrderedDict([('somecol', 'Some Value'),
                                       ('anothercol', 12345)
                                       ]),
            'type': 'Feature'}

        c.write(record)