Skip to content

Instantly share code, notes, and snippets.

@sgillies

sgillies/about.md

Last active Apr 24, 2021
Embed
What would you like to do?
Rasterio example
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
# python_powered.py
#
# Turns the 15 features in the Python logo
#
# http://www.python.org/community/logos/python-powered-h-140x182.png
#
# into a GeoJSON feature collection.
import json
import subprocess
import rasterio
from rasterio import features
with rasterio.drivers():
with rasterio.open('python-powered-h-140x182.png') as src:
blue = src.read_band(3)
# Set every non-background pixel to 0 and then mask out the
# white background.
blue[blue < 255] = 0
mask = blue != 255
# transform to world coordinates so that we can map it
transform = [-35.0, 0.5, 0.0, 40.5, 0.0, -0.5]
# Get shapes from the positive part of the image.
shapes = list(features.shapes(blue, mask=mask, transform=transform))
# Burn those shapes into a new image.
image = features.rasterize(
((g, 255) for g, v in shapes),
out_shape=blue.shape, fill=0, transform=transform )
# Write the new image out.
kwargs = src.meta
kwargs['count'] = 1
kwargs['driver'] = 'GTiff'
with rasterio.open('result.tif', 'w', **kwargs) as dst:
dst.write_band(1, image)
subprocess.call(['open', 'result.tif'])
# Write the shapes out to GeoJSON.
results = ({
'type': 'Feature',
'properties': {'raster_val': v},
'geometry': s }
for i, (s, v)
in enumerate(shapes) )
collection = {
'type': 'FeatureCollection',
'features': list(results) }
with open('python-powered.json', 'w') as dst:
json.dump(collection, dst)
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
@jseppi

This comment has been minimized.

Copy link

@jseppi jseppi commented Jan 31, 2014

Did you use a world file or some other projection source for your source image? For me, GDAL reports ERROR 6: No translation an empty SRS to PROJ.4 format is known. and I get a pretty whacky geojson representation.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment