1

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.

2

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)

Result

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)

Multiple

  • 1
    Brilliant, that did exactly what I was hoping, thanks for the quick and informative response. and I assume if I had a variable that was instead polygons I could iterate over that variable, create a record for each polygon, and do c.write(record) to put multiple rows into that shapefile? – clifgray Aug 9 '18 at 18:51
  • @clifgray, yes see edit – user2856 Aug 9 '18 at 23:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.