Hi, Jeff! I should have walked down the hall to ask you! :)

No go on closing first. fiona.open returns an unconsumed iterator so that it can handle huge featuresets. If you close it before iterating over it, it raises a ValueError "I/O operation on closed collection." The data is not actually "open in two places." It's only open on the fiona collection which is normally used within a context manager where it stays open until the collection iterator is exhausted. Look at the examples in the fiona docs.

I'll look at mapnick. But what I'm really trying to do is put together an ipython notebook showing full vector to raster trip using fiona and rasterio. The docstring for rasterio.features.rasterize advertizes this:

":param shapes: an iterator over Fiona style geometry objects (with a default value of 255) or an iterator over (geometry, value) pairs. Values must be unsigned integer type (uint8)."

And I think right there I see my problem. If "values" in the statement, "Values must be unsigned integer type (uint8)," means the feature properties, then I have a problem because:

>>> shp.schema['properties']
OrderedDict([(u'vmax', 'float'), (u'vmin', 'float')])
>>>

Looks like I need to walk my shapes and change those ranges of floats into uint8 categories. Now I understand why Sean told me yesterday that, "Rasterio isn't map-making software so this feature will always be cartographically limited..." If the data resolution is limited to 8 bits, then yeah, cartographically limited. I should use mapnick. The rasterio use case seems to be about decorating pre-existing rasters with shapes.

Just for a sanity check, I consumed the collection returned by fiona by wrapping it in list(), closed the shapefile, and then fed the list of features to rasterize, without first changing the OrderedDict to a uint8:

import fiona
from rasterio import features
collection = fiona.open("myfeatures.shp")
shapes = list(collection)
image = features.rasterize(shapes, out_shape=(1000,1000))

And that resulted in ValueError "Unsupported geometry type Feature" because my shapefile represents a collection of features where each feature is a single range for a value (i.e., four features, one for each of the ranges 0.0 to 1.0, 1.0 to 2.0, 2,0 to 3.0, and 3.0 to 4.0). I should walk the features and burn each into its own band.

But if I rasterize just one of the features with, say, shape[0], then I'm back to the seg fault, presumably because of the bad value type.

Walking the features to convert the data values to uint8 is going to be a pain because each feature is a collection of multipolygons (polygons with holes, with islands in the holes, with holes in the islands, and so on recursively). I sure hope mapnick can handle multipolygons. Even better would be if it can handle arbitrary data types other than scalars, or at least handle more data resolution than 8 bits.

Thank you for the mapnick suggestion.

--
Sincerely,

Chris Calloway, Applications Analyst
UNC Renaissance Computing Institute
100 Europa Drive, Suite 540, Chapel Hill, NC 27517
(919) 599-3530

On 10/14/2014 6:35 PM, Jeff Heard wrote:
Chris, I'd close the source dataset first.  OGR/GDAL can have serious
issues with having data open in two places, and sometimes it's
surprising what causes it.

You might also try outputting via mapnik2

On Tue, Oct 14, 2014 at 6:35 PM, Chris Calloway <c...@unc.edu
<mailto:c...@unc.edu>> wrote:

    I'm trying to programmatically generate image data from a shapefile
    with something like this:

    import fiona
    from rasterio import features
    collection = fiona.open("myfeatures.shp")
    image = features.rasterize(collection, out_shape=(1000,1000))

    And nearly immediately I get:

    "Segmentation fault: 11"

    I'm using fiona and rasterio from anaconda on OSX 10.9.5:

    $ conda list | grep rasterio
    rasterio                  0.8                      py27_0
    $ conda list | grep fiona
    fiona                     1.1.6                np19py27_0

    Has anyone gone full trip from shapefile to image data through fiona
    and rasterio.features.rasterize? If so, do you have advice? I would
    expect this is just wrapping gdal's rasterize.

    --
    Sincerely,

    Chris Calloway, Applications Analyst
    UNC Renaissance Computing Institute
    100 Europa Drive, Suite 540, Chapel Hill, NC 27517
    (919) 599-3530 <tel:%28919%29%20599-3530>

_______________________________________________
Community mailing list
Community@lists.gispython.org
http://lists.gispython.org/mailman/listinfo/community

Reply via email to