Skip to content
Browse files

big API naming cleanup

This is a big re-factor patch that tries to rationalise some of the stranger naming that has gone on. It's resulted in
some pretty chunky modules which I might have to chop down again at some point.
  • Loading branch information...
1 parent 0a5574c commit 4acd9704d5d66fe60094882eacc2bc4bea297565 @rjw57 committed Apr 25, 2012
View
8 docs-source/reference.rst
@@ -7,5 +7,11 @@ Reference
.. automodule:: foldbeam.graph
:members:
-.. automodule:: foldbeam.transform
+.. automodule:: foldbeam.raster
+ :members:
+
+.. automodule:: foldbeam.vector
+ :members:
+
+.. automodule:: foldbeam.tilestache
:members:
View
22 examples/pipeline.json
@@ -1,40 +1,32 @@
{
"nodes": {
"tilestache": {
- "type": "foldbeam.nodes:TileStacheNode",
+ "type": "foldbeam.raster:TileStacheSource",
"parameters": {
- "config_file": "base-layers-cfg.json"
+ "config_file": "../examples/base-layers-cfg.json"
}
},
"composite1": {
- "type": "foldbeam.nodes:LayerRasterNode"
+ "type": "foldbeam.raster:CompositeOver"
},
"composite2": {
- "type": "foldbeam.nodes:LayerRasterNode",
+ "type": "foldbeam.raster:CompositeOver",
"parameters": {
"top_opacity": 0.5
}
},
- "aerial": {
- "type": "foldbeam.nodes:TileStacheRasterNode"
- },
- "road": {
- "type": "foldbeam.nodes:TileStacheRasterNode"
- },
"dem": {
- "type": "foldbeam.nodes:GDALDatasetRasterNode",
+ "type": "foldbeam.raster:GDALDatasetRasterNode",
"parameters": {
"dataset": "dtm21.tif"
}
}
},
"edges": [
- [ "tilestache:yahoo_aerial", "aerial:layer" ],
- [ "tilestache:yahoo_road", "road:layer" ],
- [ "aerial:output", "composite1:bottom" ],
+ [ "tilestache:yahoo_aerial", "composite1:bottom" ],
[ "dem:output", "composite1:top" ],
[ "composite1:output", "composite2:bottom" ],
- [ "road:output", "composite2:top" ]
+ [ "tilestache:yahoo_road", "composite2:top" ]
],
"outputs": {
"output": "composite2:output"
View
4 examples/test_pipeline.py
@@ -81,7 +81,7 @@ def output_node(node, name, nodes, pads, enclosing_nodes, current_nodes):
type_color=type_colors[pad.type],
pad_name=pad_name,
pad_title=pad_title,
- align='LEFT' if pad in inputs else 'RIGHT'))
+ align='LEFT' if pad_name in node.inputs else 'RIGHT'))
pads[pad] = ('"%s":%s' % (node_name, 'pad_' + pad_name), node_name)
output.write('</TABLE>\n>\n')
@@ -143,7 +143,7 @@ def output_node(node, name, nodes, pads, enclosing_nodes, current_nodes):
output.write('}\n')
def main():
- config = json.load(open('pipeline.json'))
+ config = json.load(open(os.path.join(os.path.dirname(__file__), 'pipeline.json')))
pipeline = Pipeline(config)
pipeline_to_dot(pipeline.nodes, pipeline.outputs.values()[0], open('pipeline.dot', 'w'))
View
16 examples/test_vector_pipeline.py
@@ -15,17 +15,15 @@ def main():
config = dict(
nodes = dict(
tilestache = dict(
- type = 'foldbeam.nodes:TileStacheNode',
- parameters = dict(config_file = 'base-layers-cfg.json'),
+ type = 'foldbeam.raster:TileStacheSource',
+ parameters = dict(config_file = os.path.join(os.path.dirname(__file__), 'base-layers-cfg.json')),
),
- road = dict(type = 'foldbeam.nodes:TileStacheRasterNode'),
- aerial = dict(type = 'foldbeam.nodes:TileStacheRasterNode'),
hybrid = dict(
- type = 'foldbeam.nodes:LayerRasterNode',
+ type = 'foldbeam.raster:CompositeOver',
parameters = dict(top_opacity = 0.5),
),
composite = dict(
- type = 'foldbeam.nodes:LayerRasterNode',
+ type = 'foldbeam.raster:CompositeOver',
),
vector = dict(
type = 'foldbeam.vector:VectorRendererNode',
@@ -37,10 +35,8 @@ def main():
),
),
edges = [
- [ 'tilestache:yahoo_road', 'road:layer' ],
- [ 'tilestache:yahoo_aerial', 'aerial:layer' ],
- [ 'road:output', 'hybrid:top' ],
- [ 'aerial:output', 'hybrid:bottom' ],
+ [ 'tilestache:yahoo_road', 'hybrid:top' ],
+ [ 'tilestache:yahoo_aerial', 'hybrid:bottom' ],
[ 'hybrid:output', 'composite:bottom' ],
[ 'vector:output', 'composite:top' ],
],
View
169 src/foldbeam/core.py
@@ -211,172 +211,3 @@ def __repr__(self):
def __str__(self):
return '(%f => %f, %f => %f)' % (self.left, self.right, self.top, self.bottom)
-
-def to_rgba_unknown(array):
- """Default conversion from an array of unknown type to an RGBA array. This simply fills the array with a default
- pattern.
-
- """
-
- array = np.atleast_2d(array)
- rgba = np.empty(array.shape[:2] + (4,))
- red = np.reshape(np.arange(array.shape[1], dtype=np.float32) / array.shape[1], (1, array.shape[1]))
- green = np.reshape(np.arange(array.shape[0], dtype=np.float32) / array.shape[0], (array.shape[0], 1))
- alpha = 1
- mask = np.ma.getmask(array)
- if mask is not np.ma.nomask:
- alpha = np.where(np.any(mask, 2), 0.0, 1.0)
- rgba[:, :, 0] = np.repeat(red, array.shape[0], 0) * alpha
- rgba[:, :, 1] = np.repeat(green, array.shape[1], 1) * alpha
- rgba[:, :, 2] = 0
- rgba[:, :, 3] = alpha
- return rgba
-
-class RgbaFromBands(object):
- # Band interpretations
- RED = 'RED'
- GREEN = 'GREEN'
- BLUE = 'BLUE'
- ALPHA = 'ALPHA'
- GRAY = 'GRAY'
- NONE = 'NONE'
-
- def __init__(self, bands, is_premultiplied):
- self.bands = bands
- self.is_premultiplied = is_premultiplied
-
- def __call__(self, array):
- rgba = to_rgba_unknown(array)
- if np.any(rgba[:,:,3] != 1.0):
- mask_alpha = rgba[:,:,3]
- else:
- mask_alpha = 1.0
-
- for idx, band in enumerate(self.bands):
- scale = (band[1] if len(band) >= 2 else 1.0) * mask_alpha
- interp = band[0]
-
- if interp is RgbaFromBands.GRAY:
- rgba[:, :, :3] = np.repeat(array[:,:,idx], 3, 2) * scale
- elif interp is RgbaFromBands.RED:
- rgba[:, :, 0] = array[:,:,idx]
- rgba[:, :, 0] *= scale
- elif interp is RgbaFromBands.GREEN:
- rgba[:, :, 1] = array[:,:,idx]
- rgba[:, :, 1] *= scale
- elif interp is RgbaFromBands.BLUE:
- rgba[:, :, 2] = array[:,:,idx]
- rgba[:, :, 2] *= scale
- elif interp is RgbaFromBands.ALPHA:
- rgba[:, :, 3] = array[:,:,idx]
- rgba[:, :, 3] *= scale
-
- if not self.is_premultiplied:
- for i in xrange(3):
- rgba[:, :, i] *= rgba[:,:,3]
-
- return rgba
-
-class Raster(object):
- def __init__(self, array, envelope, to_rgba=None, can_interpolate=True, prototype=None):
- self.array = np.atleast_3d(np.float32(array))
- self.envelope = envelope
- if to_rgba is None:
- to_rgba = to_rgba_unknown
- self.to_rgba_cb = to_rgba
- self.can_interpolate = can_interpolate
-
- if prototype is not None:
- self.to_rgba_cb = prototype.to_rgba_cb
- self.can_interpolate = prototype.can_interpolate
-
- def to_rgba(self):
- return self.to_rgba_cb(self.array)
-
- def as_rgba_dataset(self):
- arr = self.to_rgba()
- if len(arr.shape) > 2:
- arr = arr.transpose((2,0,1))
- ds = gdal_array.OpenArray(np.uint8(255.0*np.clip(arr,0,1)))
- ds.SetProjection(self.envelope.spatial_reference.ExportToWkt())
- size = [ds.RasterXSize, ds.RasterYSize]
- xscale, yscale = [float(x[0])/float(x[1]) for x in zip(self.envelope.offset(), size)]
- ds.SetGeoTransform((self.envelope.left, xscale, 0, self.envelope.top, 0, yscale))
-
- return ds
-
- def as_dataset(self):
- arr = self.array
- if len(arr.shape) > 2:
- arr = arr.transpose((2,0,1))
-
- ds = gdal_array.OpenArray(arr)
-
- ds.SetProjection(self.envelope.spatial_reference.ExportToWkt())
- size = [ds.RasterXSize, ds.RasterYSize]
- xscale, yscale = [float(x[0])/float(x[1]) for x in zip(self.envelope.offset(), size)]
- ds.SetGeoTransform((self.envelope.left, xscale, 0, self.envelope.top, 0, yscale))
-
- mask = np.ma.getmask(self.array)
- mask_ds = None
- if mask is not np.ma.nomask:
- # There is some mask we need to overlay onto the dataset
- for i in xrange(1, ds.RasterCount+1):
- band = ds.GetRasterBand(i)
- band.SetNoDataValue(float('nan'))
- band.WriteArray(np.where(mask[:,:,i-1], float('nan'), self.array[:,:,i-1]))
-
- return ds
-
- def write_tiff(self, filename):
- driver = gdal.GetDriverByName('GTiff')
- driver.CreateCopy(filename, self.as_rgba_dataset())
-
- @classmethod
- def from_dataset(cls, ds, mask_band=None, **kwargs):
- ds_array = ds.ReadAsArray()
- if len(ds_array.shape) > 2:
- ds_array = ds_array.transpose((1,2,0))
- else:
- ds_array = np.atleast_3d(ds_array)
-
- if mask_band is not None:
- mask = np.atleast_3d(mask_band.ReadAsArray()) == 0
- if np.any(mask):
- ds_array = np.ma.masked_where(np.repeat(mask, ds_array.shape[2], 2), ds_array)
-
- srs = osr.SpatialReference()
- srs.ImportFromWkt(ds.GetProjection())
- envelope = dataset_envelope(ds, srs)
-
- can_interpolate = gdal.GCI_PaletteIndex not in [
- ds.GetRasterBand(i).GetColorInterpretation()
- for i in xrange(1, ds.RasterCount+1)
- ]
- if 'can_interpolate' in kwargs:
- can_interpolate = can_interpolate and kwargs['can_interpolate']
- del kwargs['can_interpolate']
-
- return Raster(ds_array, envelope, can_interpolate=can_interpolate, **kwargs)
-
-def dataset_envelope(dataset, spatial_reference):
- gt = dataset.GetGeoTransform()
- geo_trans = np.matrix([
- [ gt[1], gt[2], gt[0] ],
- [ gt[4], gt[5], gt[3] ],
- ])
-
- bounds = geo_trans * np.matrix([
- [ 0, 0, 1],
- [ 0, dataset.RasterYSize, 1 ],
- [ dataset.RasterXSize, 0, 1],
- [ dataset.RasterXSize, dataset.RasterYSize, 1 ],
- ]).transpose()
-
- min_bound = bounds.min(1)
- max_bound = bounds.max(1)
-
- return Envelope(
- bounds[0,0], bounds[0,3],
- bounds[1,0], bounds[1,3],
- spatial_reference)
View
1 src/foldbeam/editor/app.py
@@ -1,6 +1,5 @@
import logging
-from foldbeam.nodes import LayerRasterNode
from foldbeam.tilestache import TileStacheServerNode
from foldbeam.graph import node_classes, FloatType
from twisted.internet import gtk2reactor
View
3 src/foldbeam/graph.py
@@ -366,8 +366,5 @@ class NamedType(object):
def get_description(cls):
return cls.__name__
-class RasterType(NamedType):
- pass
-
class FloatType(NamedType):
pass
View
406 src/foldbeam/nodes.py
@@ -1,406 +0,0 @@
-from __future__ import print_function
-
-import logging
-import math
-
-from ModestMaps.Core import Point, Coordinate
-from osgeo import gdal
-from osgeo.osr import SpatialReference
-import numpy as np
-import TileStache
-
-from . import _gdal, core, graph, pads, transform
-from .graph import connect, ConstantNode, node
-
-#@node
-#class ToRgbaRasterNode(graph.Node):
-# def __init__(self, input_pad):
-# super(ToRgbaRasterNode, self).__init__()
-# self.add_output('output', graph.RasterType, self._render)
-# self.add_input('input', graph.RasterType)
-#
-# def _render(self, envelope, size):
-# if size is None:
-# size = map(int, envelope.size())
-#
-# raster = self.inputs.input(envelope, size)
-# if raster is None:
-# return None
-#
-# return core.Raster(raster.to_rgba(), envelope, to_rgba=lambda x: x)
-
-@node
-class LayerRasterNode(graph.Node):
- def __init__(self, top=None, bottom=None, top_opacity=None, bottom_opacity=None):
- super(LayerRasterNode, self).__init__()
- self.add_output('output', graph.RasterType, self._render)
- self.add_input('top', graph.RasterType, top)
- self.add_input('top_opacity', graph.FloatType, top_opacity)
- self.add_input('bottom', graph.RasterType, top)
- self.add_input('bottom_opacity', graph.FloatType, bottom_opacity)
-
- for input_pad in self.inputs.values():
- input_pad.damaged.connect(self._inputs_damaged)
- input_pad.connected.connect(self._inputs_damaged)
-
- def _inputs_damaged(self, boundary):
- self.outputs.output.damaged(boundary)
-
- def _render(self, envelope, size):
- opacities = [x if x is not None else 1.0 for x in [self.bottom_opacity, self.top_opacity]]
- layers = [
- self.inputs.bottom(envelope=envelope, size=size),
- self.inputs.top(envelope=envelope, size=size),
- ]
-
- output = None
- for raster, opacity in zip(layers, opacities):
- if raster is None:
- continue
- layer = raster.to_rgba()
- if layer is None:
- print('layer failed to convert')
- continue
-
- if output is None:
- output = layer
- output *= opacity
- continue
-
- one_minus_alpha = np.atleast_3d(1.0 - opacity * layer[:,:,3])
-
- output[:,:,:3] *= np.repeat(one_minus_alpha, 3, 2)
- output[:,:,:3] += opacity * layer[:,:,:3]
-
- output[:,:,3] *= one_minus_alpha[:,:,0]
- output[:,:,3] += opacity * layer[:,:,3]
-
- if output is None:
- return None
-
- return core.Raster(output, envelope, to_rgba=lambda x: x)
-
-@node
-class FileReaderNode(graph.Node):
- def __init__(self, filename=None):
- super(FileReaderNode, self).__init__()
- self.contents = None
- self.add_input('filename', str, filename)
- self.add_output('contents', gdal.Dataset, self._load)
-
- def _load(self):
- if self.contents is not None:
- return self.contents
-
- filename = self.inputs.filename
- if filename is None:
- return None
-
- self.contents = open(filename).read()
- return self.contents
-
-@node
-class GDALDatasetSourceNode(graph.Node):
- def __init__(self, filename=None):
- super(GDALDatasetSourceNode, self).__init__()
- self._dataset = None
- self._filename = None
-
- self.add_input('filename', str, filename)
- self.add_output('dataset', gdal.Dataset, lambda: self._dataset)
-
- self.inputs.filename.damaged.connect(self._filename_damaged)
- self.inputs.filename.connected.connect(self._filename_damaged)
- self._filename_damaged(None)
-
- def _filename_damaged(self, boundary, **kwargs):
- filename = self.inputs.filename()
- if filename == self._filename:
- return
-
- if filename is not None:
- logging.info('Opening GDAL dataset: ' + str(filename))
- self._dataset = gdal.Open(filename)
- logging.info('Opened GDAL dataset ' + str(self._dataset))
- else:
- logging.info('Dropping GDAL dataset ' + str(self._dataset))
- self._dataset = None
-
- self._filename = filename
-
-@node
-class GDALDatasetRasterNode(graph.Node):
- def __init__(self, dataset=None):
- super(GDALDatasetRasterNode, self).__init__()
-
- self.add_input('dataset', gdal.Dataset)
- if isinstance(dataset, basestring):
- ds_node = self.add_subnode(GDALDatasetSourceNode(dataset))
- connect(ds_node.outputs.dataset, self.inputs.dataset)
- elif dataset is not None:
- ds_node = self.add_subnode(ConstantNode(gdal.Dataset, dataset))
- connect(ds_node.outputs.dataset, self.inputs.dataset)
-
- self.add_output('output', graph.RasterType, self._render_reprojected)
-
- def _to_rgba(self, array):
- rgba = core.to_rgba_unknown(array)
- if np.any(rgba[:,:,3] != 1.0):
- mask_alpha = rgba[:,:,3]
- else:
- mask_alpha = 1.0
-
- for i in xrange(1, self.dataset.RasterCount+1):
- band = self.dataset.GetRasterBand(i)
- interp = band.GetColorInterpretation()
- data = array[:,:,i-1] * mask_alpha
-
- if interp == gdal.GCI_RedBand:
- rgba[:,:,0] = data / 255.0
- elif interp == gdal.GCI_GreenBand:
- rgba[:,:,1] = data / 255.0
- elif interp == gdal.GCI_BlueBand:
- rgba[:,:,2] = data / 255.0
- elif interp == gdal.GCI_AlphaBand:
- rgba[:,:,3] = data / 255.0
- elif interp == gdal.GCI_GrayIndex:
- rgba[:,:,:3] = np.repeat(np.atleast_3d(data / 255.0), 3, 2)
- elif interp == gdal.GCI_PaletteIndex:
- table = band.GetColorTable()
- entries = np.array([tuple(table.GetColorEntry(i)) for i in xrange(table.GetCount())])
-
- # FIXME: This assumes the palette is RGBA
- rgba = np.float32(entries[np.int32(data)]) / 255.0
- for i in xrange(4):
- rgba[:,:,i] *= mask_alpha
-
- return rgba
-
- def _render_reprojected(self, **kwargs):
- if self.dataset is None:
- return None
-
- self.spatial_reference = SpatialReference()
- self.spatial_reference.ImportFromWkt(self.dataset.GetProjection())
- self.envelope = _gdal.dataset_envelope(self.dataset, self.spatial_reference)
- self.boundary = core.boundary_from_envelope(self.envelope)
- self.is_palette = self.dataset.GetRasterBand(1).GetColorInterpretation() == gdal.GCI_PaletteIndex
-
- return pads.ReprojectingRasterFilter(
- self.spatial_reference,
- pads.TiledRasterFilter(self._render, tile_size=256))(**kwargs)
-
- def _render(self, envelope, size):
- if self.dataset is None:
- return None
-
- assert envelope.spatial_reference.IsSame(self.spatial_reference)
-
- # check if the requested area is contained within the dataset bounds
- requested_boundary = core.boundary_from_envelope(envelope)
- if not self.boundary.geometry.Intersects(requested_boundary.geometry):
- # early out if the dataset is nowhere near the requested envelope
- return None
-
- # Get the destination raster
- raster = _gdal.create_render_dataset(envelope, size, prototype_ds=self.dataset)
- ds = raster.dataset
-
- desired_srs_wkt = envelope.spatial_reference.ExportToWkt()
- gdal.ReprojectImage(
- self.dataset, ds,
- self.dataset.GetProjection(),
- desired_srs_wkt,
- gdal.GRA_NearestNeighbour if self.is_palette else gdal.GRA_Bilinear)
-
- # Create a mask raster
- mask_raster = _gdal.create_render_dataset(envelope, size, data_type=gdal.GDT_Float32)
- mask_ds = mask_raster.dataset
- band = mask_ds.GetRasterBand(1)
- band.SetColorInterpretation(gdal.GCI_Undefined)
- band.SetNoDataValue(float('nan'))
- band.Fill(float('nan'))
- gdal.ReprojectImage(
- self.dataset, mask_ds,
- self.dataset.GetProjection(),
- desired_srs_wkt,
- gdal.GRA_NearestNeighbour)
- mask_band = mask_ds.GetRasterBand(1).GetMaskBand()
-
- return core.Raster.from_dataset(ds, mask_band=mask_band, to_rgba=self._to_rgba)
-
-@node
-class TileStacheNode(graph.Node):
- def __init__(self, config_file=None):
- super(TileStacheNode, self).__init__()
- self.add_input('config_file', str, config_file)
- self.inputs.config_file.damaged.connect(self._config_updated)
- self.inputs.config_file.connected.connect(self._config_updated)
- self._config_updated()
-
- def _config_updated(self, *args, **kwargs):
- filename = self.config_file
- self.config = None
- if filename is not None:
- self.config = TileStache.parseConfigfile(filename)
- for name in sorted(self.config.layers.keys()):
- self.add_output(name, TileStache.Core.Layer, self._layer_function(name))
-
- def _layer_function(self, name):
- return lambda: self.config.layers[name]
-
-@node
-class TileStacheRasterNode(graph.Node):
- def __init__(self, layer=None):
- super(TileStacheRasterNode, self).__init__()
- self.add_input('layer', TileStache.Core.Layer, layer)
- self.add_output('output', graph.RasterType, pads.TiledRasterFilter(self._render, tile_size=256))
-
- def _update(self):
- self.preferred_srs = SpatialReference()
- self.preferred_srs.ImportFromProj4(self.layer.projection.srs)
- self.preferred_srs_wkt = self.preferred_srs.ExportToWkt()
-
- # Calculate the bounds of the zoom level 0 tile
- bounds = [
- self.layer.projection.coordinateProj(Coordinate(0,0,0)),
- self.layer.projection.coordinateProj(Coordinate(1,1,0)),
- ]
-
- self.proj_axes = (bounds[1].x - bounds[0].x, bounds[1].y - bounds[0].y)
- self.proj_size = map(abs, self.proj_axes)
- self.proj_origin = (bounds[0].x, bounds[0].y)
-
- min_x = min([p.x for p in bounds])
- max_x = max([p.x for p in bounds])
- min_y = min([p.y for p in bounds])
- max_y = max([p.y for p in bounds])
-
- self.proj_bounds = (min_x, min_y, max_x-min_x, max_y-min_y)
-
- def _zoom_for_envelope(self, envelope, size):
- # How many tiles should cover each axis
- n_tiles = map(lambda x: x/256.0, size)
-
- # Over what range?
- envelope_size = map(abs, envelope.offset())
- proj_range = min(envelope_size)
-
- # Giving a rough tile size in projection coordinate of...
- zoomed_tile_size = map(lambda x: x[0]/x[1], zip(envelope_size, n_tiles))
-
- # And convert to a zoom
- zooms = map(lambda x: math.log(x[0], 2) - math.log(x[1], 2), zip(self.proj_size, zoomed_tile_size))
-
- # Which is closer to the precise zoom we want?
- zoom = max(zooms)
- floor_diff = abs(math.pow(2, zoom) - math.pow(2, math.floor(zoom)))
- ceil_diff = abs(math.pow(2, zoom) - math.pow(2, math.ceil(zoom)))
-
- # Choose an overall zoom from these
- int_zoom = int(math.ceil(zoom)) if ceil_diff < floor_diff else int(math.floor(zoom))
- return max(0, min(18, int_zoom))
-
- def _render(self, envelope, size):
- if self.layer is None:
- return None
-
- if size is None:
- size = [int(x) for x in envelope.size()]
-
- # FIXME: Once output hints are enabled, they should be used here
- self._update()
- try:
- pref_envelope = envelope.transform_to(
- self.preferred_srs,
- min(envelope.size()) / float(max(size)))
- except core.ProjectionError:
- # Skip projection errors
- return None
-
- zoom = self._zoom_for_envelope(pref_envelope, size)
-
- # Get the destination raster
- raster = _gdal.create_render_dataset(envelope, size, 4)
- assert size == (raster.dataset.RasterXSize, raster.dataset.RasterYSize)
-
- # Convert the envelope into the preferred spatial reference
- try:
- pref_envelope = envelope.transform_to(
- self.preferred_srs,
- min(envelope.size()) / float(max(size)))
- except core.ProjectionError:
- print('Ignoring projection error and returning empty raster')
- return raster
-
- # How many tiles overall in x and y at this zoom
- n_proj_tiles = math.pow(2, zoom)
-
- # Compute the tile co-ordinates of each corner
- corners = [
- self.layer.projection.projCoordinate(Point(pref_envelope.left, pref_envelope.top)),
- self.layer.projection.projCoordinate(Point(pref_envelope.right, pref_envelope.top)),
- self.layer.projection.projCoordinate(Point(pref_envelope.left, pref_envelope.bottom)),
- self.layer.projection.projCoordinate(Point(pref_envelope.right, pref_envelope.bottom)),
- ]
- corners = [c.zoomTo(zoom) for c in corners]
-
- corner_rows = [int(math.floor(x.row)) for x in corners]
- corner_columns = [int(math.floor(x.column)) for x in corners]
-
- tile_rasters = []
-
- # Get each tile image
- png_driver = gdal.GetDriverByName('PNG')
- desired_srs_wkt = envelope.spatial_reference.ExportToWkt()
- assert png_driver is not None
- for r in xrange(min(corner_rows), max(corner_rows)+1):
- if r < 0 or r >= n_proj_tiles:
- continue
- for c in xrange(min(corner_columns), max(corner_columns)+1):
- tile_tl_point = self.layer.projection.coordinateProj(Coordinate(r,c,zoom))
- tile_br_point = self.layer.projection.coordinateProj(Coordinate(r+1,c+1,zoom))
-
- c = c % n_proj_tiles
- if c < 0:
- c += n_proj_tiles
- tile_coord = Coordinate(r, c, zoom)
-
- try:
- tile_type, png_data = TileStache.getTile(self.layer, tile_coord, 'png')
- if tile_type != 'image/png':
- print('Did not get PNG data when fetching tile %s. Skipping.' % (tile_coord))
- continue
- except IOError as e:
- print('Ignoring error fetching tile %s (%s).' % (tile_coord, e))
- continue
-
- gdal.FileFromMemBuffer('/vsimem/tmptile.png', png_data)
- tile_raster = gdal.Open('/vsimem/tmptile.png')
- tile_raster.SetProjection(self.preferred_srs_wkt)
-
- xscale = (tile_br_point.x - tile_tl_point.x) / tile_raster.RasterXSize
- yscale = (tile_br_point.y - tile_tl_point.y) / tile_raster.RasterXSize
- tile_raster.SetGeoTransform((
- tile_tl_point.x, xscale, 0.0,
- tile_tl_point.y, 0.0, yscale,
- ))
-
- tile_rasters.append(core.Raster.from_dataset(tile_raster))
-
- tile_raster = None
- gdal.Unlink('/vsimem/tmptile.png')
-
- output = core.Raster(
- np.ma.masked_all((size[1], size[0], 4), dtype=np.float32),
- envelope,
- to_rgba=core.RgbaFromBands(
- (
- (core.RgbaFromBands.RED, 1.0/255.0),
- (core.RgbaFromBands.GREEN, 1.0/255.0),
- (core.RgbaFromBands.BLUE, 1.0/255.0),
- (core.RgbaFromBands.ALPHA, 1.0/255.0),
- ), False)
- )
- transform.reproject_rasters(output, tile_rasters)
- return output
View
124 src/foldbeam/pads.py
@@ -1,124 +0,0 @@
-from . import _gdal, core, graph
-from .graph import InputPad, OutputPad
-import numpy as np
-from osgeo import osr, gdal
-
-class TiledRasterFilter(object):
- def __init__(self, render_cb, tile_size=None):
- self.render_cb = render_cb
- self.tile_size = tile_size
-
- def __call__(self, envelope=None, size=None, **kwargs):
- if envelope is None:
- return None
-
- if size is None:
- size = map(int, envelope.size())
-
- if self.tile_size is None:
- raster = self.render_cb(envelope=envelope, size=size, **kwargs)
- if raster is None or len(raster) == 0 or raster[0] is None:
- return None
- return raster[0]
-
- tiles = []
- tile_offsets = []
- xscale, yscale = [x[0] / float(x[1]) for x in zip(envelope.offset(), size)]
- for x in xrange(0, size[0], self.tile_size):
- width = min(x + self.tile_size, size[0]) - x
- for y in xrange(0, size[1], self.tile_size):
- height = min(y + self.tile_size, size[1]) - y
- tile_envelope = core.Envelope(
- envelope.left + xscale * x,
- envelope.left + xscale * (x + width),
- envelope.top + yscale * y,
- envelope.top + yscale * (y + height),
- envelope.spatial_reference,
- )
- tiles.append(self.render_cb(envelope=tile_envelope, size=(width, height), **kwargs))
- tile_offsets.append((x,y))
-
- results = [x for x in zip(tile_offsets, tiles) if x[1] is not None]
- if len(results) == 0:
- return None
-
- # FIXME: This assumes that to_rgba_cb is the same for each tile
- prototype = results[0][1]
- depth = prototype.array.shape[2]
-
- shape = (size[1], size[0])
- mask = np.ones(shape + (depth,), dtype=np.bool)
- data = np.zeros(shape + (depth,), dtype=np.float32)
- for pos, tile in results:
- x, y = pos
- h, w = tile.array.shape[:2]
-
- tile_mask = np.ma.getmask(tile.array)
- if tile_mask is np.ma.nomask:
- tile_mask = False
-
- mask[y:(y+h), x:(x+w), :] = tile_mask
- data[y:(y+h), x:(x+w), :] = tile.array
-
- if np.any(mask):
- data = np.ma.array(data, mask=mask)
-
- return core.Raster(data, envelope, prototype=prototype)
-
-class ReprojectingRasterFilter(object):
- def __init__(self, native_spatial_reference, render_cb):
- self.native_spatial_reference = native_spatial_reference
- self.native_spatial_reference_wkt = native_spatial_reference.ExportToWkt()
- self.render_cb = render_cb
-
- def __call__(self, envelope=None, size=None, **kwargs):
- if envelope is None:
- return None
-
- if size is None:
- size = map(int, envelope.size())
-
- if envelope.spatial_reference.IsSame(self.native_spatial_reference):
- return self.render_cb(envelope=envelope, size=size, **kwargs)
-
- # We need to reproject this data. Convert the envelope into the native spatial reference
- try:
- native_envelope = envelope.transform_to(
- self.native_spatial_reference,
- min(envelope.size()) / float(max(size)))
- except core.ProjectionError:
- # If we fail to reproject, return a null tile
- return None
-
- # Get the native tile
- raster = self.render_cb(envelope=native_envelope, size=size, **kwargs)
- if raster is None:
- return None
-
- raster_ds = raster.as_dataset()
-
- # Get the destination raster
- ds = _gdal.create_render_dataset(envelope, size, 4).dataset
-
- desired_srs_wkt = envelope.spatial_reference.ExportToWkt()
- gdal.ReprojectImage(
- raster_ds, ds,
- self.native_spatial_reference_wkt,
- desired_srs_wkt,
- gdal.GRA_Bilinear if raster.can_interpolate else gdal.GRA_NearestNeighbour)
-
- # Create a mask raster
- mask_raster = _gdal.create_render_dataset(envelope, size, 1, data_type=gdal.GDT_Float32)
- mask_ds = mask_raster.dataset
- band = mask_ds.GetRasterBand(1)
- band.SetColorInterpretation(gdal.GCI_Undefined)
- band.SetNoDataValue(float('nan'))
- band.Fill(float('nan'))
- gdal.ReprojectImage(
- raster_ds, mask_ds,
- self.native_spatial_reference_wkt,
- desired_srs_wkt,
- gdal.GRA_NearestNeighbour)
- mask_band = mask_ds.GetRasterBand(1).GetMaskBand()
-
- return core.Raster.from_dataset(ds, mask_band=mask_band, prototype=raster)
View
5 src/foldbeam/pipeline.py
@@ -1,7 +1,10 @@
from __future__ import print_function
-from .graph import connect, Node
import traceback
+from foldbeam.graph import connect, Node
+import foldbeam.vector
+import foldbeam.tilestache
+
class Pipeline(Node):
def __init__(self, configuration):
super(Pipeline, self).__init__()
View
848 src/foldbeam/raster.py
@@ -0,0 +1,848 @@
+"""
+Handling raster data
+====================
+
+"""
+from __future__ import print_function
+
+import logging
+import math
+import os
+
+from ModestMaps.Core import Point, Coordinate
+import numpy as np
+from osgeo import gdal, gdal_array, osr
+import pyproj
+import TileStache
+
+from foldbeam import _gdal, core, graph
+
+class Raster(object):
+ """An array of pixels associated with a rectangle in some geographic projection. A raster is a wrapper around a
+ :py:class:`numpy.array` instance which can be accessed directly via the :py:attr:`array` attribute. A raster has an
+ associated envelope and the left, right, top and bottom edges of the array are mapped to the left, right, top and
+ bottom of the envelope.
+
+ :param array: the initial contents of the raster
+ :type array: :py:class:`numpy.array` or similar
+ :param envelope: the bounding box and associated projection which encloses this raster
+ :type envelope: :py:class:`core.Envelope`
+ :param to_rgba: called with :py:attr:`array` to convert this raster into a 4-channel RGBA raster
+ :type to_rgba: callable
+ :param can_interpolate: flag indicating whether interpolation of values in this raster is valid
+ :param prototype: copy *can_interpolate* and *to_rgba* from a prototype :py:class:`Raster`
+ :type prototype: :py:class:`Raster`
+
+ """
+
+ def __init__(self, array, envelope, to_rgba=None, can_interpolate=True, prototype=None):
+ self.array = np.atleast_3d(np.float32(array))
+ self.envelope = envelope
+ if to_rgba is None:
+ to_rgba = to_rgba_unknown
+ self.to_rgba_cb = to_rgba
+ self.can_interpolate = can_interpolate
+
+ if prototype is not None:
+ self.to_rgba_cb = prototype.to_rgba_cb
+ self.can_interpolate = prototype.can_interpolate
+
+ def to_rgba(self):
+ """Obtain an array with shape (height, width, 4) where the depth are respectively the red, green, blue and alpha
+ values for displaying this raster on the interval [0,1).
+
+ """
+ return self.to_rgba_cb(self.array)
+
+ def as_rgba_dataset(self):
+ """Obtain a :py:class:`gdal.Dataset` with red, green, blue and alpha bands representing the array returned from
+ :py:meth:`to_rgba`.
+
+ """
+ arr = self.to_rgba()
+ if len(arr.shape) > 2:
+ arr = arr.transpose((2,0,1))
+ ds = gdal_array.OpenArray(np.uint8(255.0*np.clip(arr,0,1)))
+ ds.SetProjection(self.envelope.spatial_reference.ExportToWkt())
+ size = [ds.RasterXSize, ds.RasterYSize]
+ xscale, yscale = [float(x[0])/float(x[1]) for x in zip(self.envelope.offset(), size)]
+ ds.SetGeoTransform((self.envelope.left, xscale, 0, self.envelope.top, 0, yscale))
+
+ return ds
+
+ def as_dataset(self):
+ """Obtain a version of this raster as a :py:class:`gdal.Dataset`."""
+
+ arr = self.array
+ if len(arr.shape) > 2:
+ arr = arr.transpose((2,0,1))
+
+ ds = gdal_array.OpenArray(arr)
+
+ ds.SetProjection(self.envelope.spatial_reference.ExportToWkt())
+ size = [ds.RasterXSize, ds.RasterYSize]
+ xscale, yscale = [float(x[0])/float(x[1]) for x in zip(self.envelope.offset(), size)]
+ ds.SetGeoTransform((self.envelope.left, xscale, 0, self.envelope.top, 0, yscale))
+
+ mask = np.ma.getmask(self.array)
+ mask_ds = None
+ if mask is not np.ma.nomask:
+ # There is some mask we need to overlay onto the dataset
+ for i in xrange(1, ds.RasterCount+1):
+ band = ds.GetRasterBand(i)
+ band.SetNoDataValue(float('nan'))
+ band.WriteArray(np.where(mask[:,:,i-1], float('nan'), self.array[:,:,i-1]))
+
+ return ds
+
+ def write_tiff(self, filename):
+ """Write this raster to a GeoTIFF."""
+ driver = gdal.GetDriverByName('GTiff')
+ driver.CreateCopy(filename, self.as_rgba_dataset())
+
+ @classmethod
+ def from_dataset(cls, ds, mask_band=None, **kwargs):
+ """Create a :py:class:`Raster` from a :py:class:`gdal.Dataset`.
+
+ Extra keyword arguments are passed to the :py:class:`Raster` constructor.
+
+ """
+ ds_array = ds.ReadAsArray()
+ if len(ds_array.shape) > 2:
+ ds_array = ds_array.transpose((1,2,0))
+ else:
+ ds_array = np.atleast_3d(ds_array)
+
+ if mask_band is not None:
+ mask = np.atleast_3d(mask_band.ReadAsArray()) == 0
+ if np.any(mask):
+ ds_array = np.ma.masked_where(np.repeat(mask, ds_array.shape[2], 2), ds_array)
+
+ srs = osr.SpatialReference()
+ srs.ImportFromWkt(ds.GetProjection())
+ envelope = dataset_envelope(ds, srs)
+
+ can_interpolate = gdal.GCI_PaletteIndex not in [
+ ds.GetRasterBand(i).GetColorInterpretation()
+ for i in xrange(1, ds.RasterCount+1)
+ ]
+ if 'can_interpolate' in kwargs:
+ can_interpolate = can_interpolate and kwargs['can_interpolate']
+ del kwargs['can_interpolate']
+
+ return Raster(ds_array, envelope, can_interpolate=can_interpolate, **kwargs)
+
+def to_rgba_unknown(array):
+ """Default conversion from an array of unknown type to an RGBA array. This simply fills the array with a default
+ pattern.
+
+ """
+
+ array = np.atleast_2d(array)
+ rgba = np.empty(array.shape[:2] + (4,))
+ red = np.reshape(np.arange(array.shape[1], dtype=np.float32) / array.shape[1], (1, array.shape[1]))
+ green = np.reshape(np.arange(array.shape[0], dtype=np.float32) / array.shape[0], (array.shape[0], 1))
+ alpha = 1
+ mask = np.ma.getmask(array)
+ if mask is not np.ma.nomask:
+ alpha = np.where(np.any(mask, 2), 0.0, 1.0)
+ rgba[:, :, 0] = np.repeat(red, array.shape[0], 0) * alpha
+ rgba[:, :, 1] = np.repeat(green, array.shape[1], 1) * alpha
+ rgba[:, :, 2] = 0
+ rgba[:, :, 3] = alpha
+ return rgba
+
+class RgbaFromBands(object):
+ # Band interpretations
+ RED = 'RED'
+ GREEN = 'GREEN'
+ BLUE = 'BLUE'
+ ALPHA = 'ALPHA'
+ GRAY = 'GRAY'
+ NONE = 'NONE'
+
+ def __init__(self, bands, is_premultiplied):
+ self.bands = bands
+ self.is_premultiplied = is_premultiplied
+
+ def __call__(self, array):
+ rgba = to_rgba_unknown(array)
+ if np.any(rgba[:,:,3] != 1.0):
+ mask_alpha = rgba[:,:,3]
+ else:
+ mask_alpha = 1.0
+
+ for idx, band in enumerate(self.bands):
+ scale = (band[1] if len(band) >= 2 else 1.0) * mask_alpha
+ interp = band[0]
+
+ if interp is RgbaFromBands.GRAY:
+ rgba[:, :, :3] = np.repeat(array[:,:,idx], 3, 2) * scale
+ elif interp is RgbaFromBands.RED:
+ rgba[:, :, 0] = array[:,:,idx]
+ rgba[:, :, 0] *= scale
+ elif interp is RgbaFromBands.GREEN:
+ rgba[:, :, 1] = array[:,:,idx]
+ rgba[:, :, 1] *= scale
+ elif interp is RgbaFromBands.BLUE:
+ rgba[:, :, 2] = array[:,:,idx]
+ rgba[:, :, 2] *= scale
+ elif interp is RgbaFromBands.ALPHA:
+ rgba[:, :, 3] = array[:,:,idx]
+ rgba[:, :, 3] *= scale
+
+ if not self.is_premultiplied:
+ for i in xrange(3):
+ rgba[:, :, i] *= rgba[:,:,3]
+
+ return rgba
+
+def dataset_envelope(dataset, spatial_reference):
+ gt = dataset.GetGeoTransform()
+ geo_trans = np.matrix([
+ [ gt[1], gt[2], gt[0] ],
+ [ gt[4], gt[5], gt[3] ],
+ ])
+
+ bounds = geo_trans * np.matrix([
+ [ 0, 0, 1],
+ [ 0, dataset.RasterYSize, 1 ],
+ [ dataset.RasterXSize, 0, 1],
+ [ dataset.RasterXSize, dataset.RasterYSize, 1 ],
+ ]).transpose()
+
+ min_bound = bounds.min(1)
+ max_bound = bounds.max(1)
+
+ return core.Envelope(
+ bounds[0,0], bounds[0,3],
+ bounds[1,0], bounds[1,3],
+ spatial_reference
+ )
+
+def compute_warp(dst_envelope, dst_size, src_spatial_reference):
+ """Compute the projection co-ordinates in a source spatial projection for pixels in a destination image.
+
+ :param dst_envelope: the envelope of the destination raster
+ :type dst_envelope: :py:class:`core.Envelope`
+ :param dst_size: the width and height of the destination raster
+ :type dst_size: pair of integer
+ :param src_spatial_reference: co-ordinate system of the source raster
+ :type src_spatial_reference: :py:class:`osgeo.osr.SpatialReference`
+ :rtype: pair of arrays
+
+ Compute the projection co-ordinates of each pixel in a destination raster when projected into the source spatial
+ reference. This is necessary to re-project images correctly.
+
+ The return value is a pair of array-like objects giving the x and y co-ordinates. The returned arrays are the width
+ and height specified in dst_size.
+
+ """
+
+ # Compute x- and y-co-ordinates in dst from the envelope
+ dst_cols, dst_rows = dst_size
+ dst_x = np.repeat(
+ np.linspace(dst_envelope.left, dst_envelope.right, dst_cols).reshape(1, dst_cols),
+ dst_rows, 0)
+ dst_y = np.repeat(
+ np.linspace(dst_envelope.top, dst_envelope.bottom, dst_rows).reshape(dst_rows, 1),
+ dst_cols, 1)
+
+ assert dst_x.shape == (dst_rows, dst_cols)
+ assert dst_y.shape == (dst_rows, dst_cols)
+
+ # Special case when source and destination have same srs
+ if dst_envelope.spatial_reference.IsSame(src_spatial_reference):
+ return dst_x, dst_y
+
+ dst_proj = pyproj.Proj(dst_envelope.spatial_reference.ExportToProj4())
+ src_proj = pyproj.Proj(src_spatial_reference.ExportToProj4())
+
+ src_x, src_y = pyproj.transform(dst_proj, src_proj, dst_x, dst_y)
+
+ return src_x, src_y
+
+def proj_to_pixels(x, y, envelope, size):
+ """Compute pixel co-ordinates of projection co-ordinates x and y in image with envelope and size specified."""
+
+ off = envelope.offset()
+ x = (x - envelope.left) * (float(size[0]-1) / off[0])
+ y = (y - envelope.top) * (float(size[1]-1) / off[1])
+ return x,y
+
+def pixels_to_proj(x, y, envelope, size):
+ """Compute projection co-ordinates of pixel co-ordinates x and y in image with envelope and size specified."""
+
+ off = envelope.offset()
+ x = x * (off[0] / float(size[0]-1)) + envelope.left
+ y = y * (off[1] / float(size[1]-1)) + envelope.right
+ return x,y
+
+def sample_raster(src_pixel_x, src_pixel_y, src_raster):
+ """Sample pixels from a raster.
+
+ :param src_pixel_x: array of pixel x-co-ordinates
+ :param src_pixel_y: array of pixel y-co-ordinates
+ :param src_raster: source raster
+ :type src_raster: :py:class:`Raster`
+ :rtype: :py:class:`numpy.array` or :py:class:`numpy.ma.array`
+
+ The returned array has the same shape as src_raster with each value corresponding to the sampled value. The arrays
+ src_pixel_x and src_pixel_y must have the same width and height as src_raster. If pixel co-ordinates outside of the
+ range of those present in src_raster are specified, the values will be masked out of the resulting array and a numpy
+ masked array will be returned.
+
+ """
+
+ px = np.int32(np.round(src_pixel_x))
+ py = np.int32(np.round(src_pixel_y))
+
+ x_valid = np.logical_and(px >= 0, px < src_raster.array.shape[1])
+ y_valid = np.logical_and(py >= 0, py < src_raster.array.shape[0])
+ valid = np.logical_and(x_valid, y_valid)
+
+ valid_flag = valid.ravel()
+ valid_src_indices = np.ravel_multi_index(
+ (py.ravel()[valid_flag], px.ravel()[valid_flag]),
+ src_raster.array.shape[:2]
+ )
+
+ dst_shape = np.atleast_2d(src_pixel_x).shape[:2] + (src_raster.array.shape[2],)
+ dst = np.empty(shape=dst_shape, dtype=np.float32)
+
+ for i in xrange(dst_shape[2]):
+ dst[valid, i] = src_raster.array[:,:,i].ravel()[valid_src_indices]
+
+ dst = np.ma.fix_invalid(
+ dst,
+ fill_value=np.nan,
+ mask=np.repeat(np.atleast_3d(np.logical_not(valid)), dst_shape[2], 2)
+ )
+
+ return dst
+
+def reproject_raster(dst, src):
+ """Re-project a source raster into a destination raster.
+
+ See :py:func:`reproject_rasters`.
+
+ """
+ reproject_rasters(dst, (src,))
+
+def _py_reproject_rasters(dst, srcs):
+ if len(srcs) == 0:
+ return
+
+ src_srs = srcs[0].envelope.spatial_reference
+ assert all([x.envelope.spatial_reference.IsSame(src_srs) for x in srcs[1:]])
+
+ dst_height, dst_width = dst.array.shape[:2]
+ src_x, src_y = compute_warp(dst.envelope, (dst_width, dst_height), src_srs)
+ for src in srcs:
+ src_height, src_width = src.array.shape[:2]
+
+ # Compute warp
+ src_px, src_py = proj_to_pixels(src_x, src_y, src.envelope, (src_width, src_height))
+ sample = sample_raster(src_px, src_py, src)
+
+ bands = min(sample.shape[2], dst.array.shape[2])
+ dst.array[:,:,:bands] = np.ma.filled(sample[:,:,:bands], dst.array[:,:,:bands])
+
+def _gdal_reproject_rasters(dst, srcs):
+ if len(srcs) == 0:
+ return
+
+ src_srs = srcs[0].envelope.spatial_reference
+ assert all([x.envelope.spatial_reference.IsSame(src_srs) for x in srcs[1:]])
+
+ mem_driver = gdal.GetDriverByName('MEM')
+ dst_ds = mem_driver.CreateCopy('', dst.as_dataset())
+ src_wkt = src_srs.ExportToWkt()
+ dst_wkt = dst.envelope.spatial_reference.ExportToWkt()
+ for src in srcs:
+ gdal.ReprojectImage(
+ src.as_dataset(), dst_ds,
+ src_wkt, dst_wkt,
+ gdal.GRA_Bilinear if src.can_interpolate else gdal.GRA_NearestNeighbour)
+
+ dst.array = Raster.from_dataset(dst_ds).array
+
+if 'FOLDBEAM_REPROJECTION_PROVIDER' in os.environ:
+ _provider = os.environ['FOLDBEAM_REPROJECTION_PROVIDER']
+ if _provider == 'GDAL':
+ reproject_rasters = _gdal_reproject_rasters
+ elif _provider == 'PROJ':
+ reproject_rasters = _py_reproject_rasters
+ else:
+ raise ValueError('Unknown projection provider: ' + _provider)
+else:
+ reproject_rasters = _gdal_reproject_rasters
+
+@graph.node
+class CompositeOver(graph.Node):
+ def __init__(self, top=None, bottom=None, top_opacity=None, bottom_opacity=None):
+ super(CompositeOver, self).__init__()
+ self.add_output('output', Raster, self._render)
+ self.add_input('top', Raster, top)
+ self.add_input('top_opacity', graph.FloatType, top_opacity)
+ self.add_input('bottom', Raster, top)
+ self.add_input('bottom_opacity', graph.FloatType, bottom_opacity)
+
+ for input_pad in self.inputs.values():
+ input_pad.damaged.connect(self._inputs_damaged)
+ input_pad.connected.connect(self._inputs_damaged)
+
+ def _inputs_damaged(self, boundary):
+ self.outputs.output.damaged(boundary)
+
+ def _render(self, envelope, size):
+ opacities = [x if x is not None else 1.0 for x in [self.bottom_opacity, self.top_opacity]]
+ layers = [
+ self.inputs.bottom(envelope=envelope, size=size),
+ self.inputs.top(envelope=envelope, size=size),
+ ]
+
+ output = None
+ for raster, opacity in zip(layers, opacities):
+ if raster is None:
+ continue
+ layer = raster.to_rgba()
+ if layer is None:
+ print('layer failed to convert')
+ continue
+
+ if output is None:
+ output = layer
+ output *= opacity
+ continue
+
+ one_minus_alpha = np.atleast_3d(1.0 - opacity * layer[:,:,3])
+
+ output[:,:,:3] *= np.repeat(one_minus_alpha, 3, 2)
+ output[:,:,:3] += opacity * layer[:,:,:3]
+
+ output[:,:,3] *= one_minus_alpha[:,:,0]
+ output[:,:,3] += opacity * layer[:,:,3]
+
+ if output is None:
+ return None
+
+ return Raster(output, envelope, to_rgba=lambda x: x)
+
+@graph.node
+class GDALDatasetSourceNode(graph.Node):
+ def __init__(self, filename=None):
+ super(GDALDatasetSourceNode, self).__init__()
+ self._dataset = None
+ self._filename = None
+
+ self.add_input('filename', str, filename)
+ self.add_output('dataset', gdal.Dataset, lambda: self._dataset)
+
+ self.inputs.filename.damaged.connect(self._filename_damaged)
+ self.inputs.filename.connected.connect(self._filename_damaged)
+ self._filename_damaged(None)
+
+ def _filename_damaged(self, boundary, **kwargs):
+ filename = self.inputs.filename()
+ if filename == self._filename:
+ return
+
+ if filename is not None:
+ logging.info('Opening GDAL dataset: ' + str(filename))
+ self._dataset = gdal.Open(filename)
+ logging.info('Opened GDAL dataset ' + str(self._dataset))
+ else:
+ logging.info('Dropping GDAL dataset ' + str(self._dataset))
+ self._dataset = None
+
+ self._filename = filename
+
+@graph.node
+class GDALDatasetRasterNode(graph.Node):
+ def __init__(self, dataset=None):
+ super(GDALDatasetRasterNode, self).__init__()
+
+ self.add_input('dataset', gdal.Dataset)
+ if isinstance(dataset, basestring):
+ ds_node = self.add_subnode(GDALDatasetSourceNode(dataset))
+ graph.connect(ds_node.outputs.dataset, self.inputs.dataset)
+ elif dataset is not None:
+ ds_node = self.add_subnode(graph.ConstantNode(gdal.Dataset, dataset))
+ graph.connect(ds_node.outputs.dataset, self.inputs.dataset)
+
+ self.add_output('output', Raster, self._render_reprojected)
+
+ def _to_rgba(self, array):
+ rgba = to_rgba_unknown(array)
+ if np.any(rgba[:,:,3] != 1.0):
+ mask_alpha = rgba[:,:,3]
+ else:
+ mask_alpha = 1.0
+
+ for i in xrange(1, self.dataset.RasterCount+1):
+ band = self.dataset.GetRasterBand(i)
+ interp = band.GetColorInterpretation()
+ data = array[:,:,i-1] * mask_alpha
+
+ if interp == gdal.GCI_RedBand:
+ rgba[:,:,0] = data / 255.0
+ elif interp == gdal.GCI_GreenBand:
+ rgba[:,:,1] = data / 255.0
+ elif interp == gdal.GCI_BlueBand:
+ rgba[:,:,2] = data / 255.0
+ elif interp == gdal.GCI_AlphaBand:
+ rgba[:,:,3] = data / 255.0
+ elif interp == gdal.GCI_GrayIndex:
+ rgba[:,:,:3] = np.repeat(np.atleast_3d(data / 255.0), 3, 2)
+ elif interp == gdal.GCI_PaletteIndex:
+ table = band.GetColorTable()
+ entries = np.array([tuple(table.GetColorEntry(i)) for i in xrange(table.GetCount())])
+
+ # FIXME: This assumes the palette is RGBA
+ rgba = np.float32(entries[np.int32(data)]) / 255.0
+ for i in xrange(4):
+ rgba[:,:,i] *= mask_alpha
+
+ return rgba
+
+ def _render_reprojected(self, **kwargs):
+ if self.dataset is None:
+ return None
+
+ self.spatial_reference = osr.SpatialReference()
+ self.spatial_reference.ImportFromWkt(self.dataset.GetProjection())
+ self.envelope = dataset_envelope(self.dataset, self.spatial_reference)
+ self.boundary = core.boundary_from_envelope(self.envelope)
+ self.is_palette = self.dataset.GetRasterBand(1).GetColorInterpretation() == gdal.GCI_PaletteIndex
+
+ return ReprojectingRasterFilter(
+ self.spatial_reference,
+ TiledRasterFilter(self._render, tile_size=256))(**kwargs)
+
+ def _render(self, envelope, size):
+ if self.dataset is None:
+ return None
+
+ assert envelope.spatial_reference.IsSame(self.spatial_reference)
+
+ # check if the requested area is contained within the dataset bounds
+ requested_boundary = core.boundary_from_envelope(envelope)
+ if not self.boundary.geometry.Intersects(requested_boundary.geometry):
+ # early out if the dataset is nowhere near the requested envelope
+ return None
+
+ # Get the destination raster
+ raster = _gdal.create_render_dataset(envelope, size, prototype_ds=self.dataset)
+ ds = raster.dataset
+
+ desired_srs_wkt = envelope.spatial_reference.ExportToWkt()
+ gdal.ReprojectImage(
+ self.dataset, ds,
+ self.dataset.GetProjection(),
+ desired_srs_wkt,
+ gdal.GRA_NearestNeighbour if self.is_palette else gdal.GRA_Bilinear)
+
+ # Create a mask raster
+ mask_raster = _gdal.create_render_dataset(envelope, size, data_type=gdal.GDT_Float32)
+ mask_ds = mask_raster.dataset
+ band = mask_ds.GetRasterBand(1)
+ band.SetColorInterpretation(gdal.GCI_Undefined)
+ band.SetNoDataValue(float('nan'))
+ band.Fill(float('nan'))
+ gdal.ReprojectImage(
+ self.dataset, mask_ds,
+ self.dataset.GetProjection(),
+ desired_srs_wkt,
+ gdal.GRA_NearestNeighbour)
+ mask_band = mask_ds.GetRasterBand(1).GetMaskBand()
+
+ return Raster.from_dataset(ds, mask_band=mask_band, to_rgba=self._to_rgba)
+
+@graph.node
+class TileStacheSource(graph.Node):
+ def __init__(self, config_file=None, config=None):
+ super(TileStacheSource, self).__init__()
+ self.add_input('config_file', str, config_file)
+ self._config = None
+ self.inputs.config_file.damaged.connect(self._config_updated)
+ self.inputs.config_file.connected.connect(self._config_updated)
+ self._config_updated()
+
+ if config is not None:
+ self.config = config
+
+ @property
+ def config(self):
+ return self._config
+
+ @config.setter
+ def config(self, config):
+ self._config = config
+ self._update_pads_from_config()
+
+ def _config_updated(self, *args, **kwargs):
+ filename = self.config_file
+ self.config = None
+ if filename is not None:
+ self.config = TileStache.parseConfigfile(filename)
+
+ def _update_pads_from_config(self):
+ if self.config is None:
+ return
+ for name in sorted(self.config.layers.keys()):
+ self.add_output(name, Raster, self._layer_function(name))
+
+ def _layer_function(self, name):
+ return lambda **kwargs: self._render(layer=self.config.layers[name], **kwargs)
+
+ def _zoom_for_envelope(self, layer, envelope, size):
+ # How many tiles should cover each axis
+ n_tiles = map(lambda x: x/256.0, size)
+
+ # Over what range?
+ envelope_size = map(abs, envelope.offset())
+ proj_range = min(envelope_size)
+
+ # Calculate the bounds of the zoom level 0 tile
+ bounds = [
+ layer.projection.coordinateProj(Coordinate(0,0,0)),
+ layer.projection.coordinateProj(Coordinate(1,1,0)),
+ ]
+ proj_axes = (bounds[1].x - bounds[0].x, bounds[1].y - bounds[0].y)
+ proj_size = map(abs, proj_axes)
+
+ # Giving a rough tile size in projection coordinate of...
+ zoomed_tile_size = map(lambda x: x[0]/x[1], zip(envelope_size, n_tiles))
+
+ # And convert to a zoom
+ zooms = map(lambda x: math.log(x[0], 2) - math.log(x[1], 2), zip(proj_size, zoomed_tile_size))
+
+ # Which is closer to the precise zoom we want?
+ zoom = max(zooms)
+ floor_diff = abs(math.pow(2, zoom) - math.pow(2, math.floor(zoom)))
+ ceil_diff = abs(math.pow(2, zoom) - math.pow(2, math.ceil(zoom)))
+
+ # Choose an overall zoom from these
+ int_zoom = int(math.ceil(zoom)) if ceil_diff < floor_diff else int(math.floor(zoom))
+ return max(0, min(18, int_zoom))
+
+ def _render(self, layer=None, envelope=None, size=None):
+ if layer is None or envelope is None:
+ return None
+
+ if size is None:
+ size = [int(x) for x in envelope.size()]
+
+ preferred_srs = osr.SpatialReference()
+ preferred_srs.ImportFromProj4(layer.projection.srs)
+ preferred_srs_wkt = preferred_srs.ExportToWkt()
+
+ # FIXME: Once output hints are enabled, they should be used here
+ try:
+ pref_envelope = envelope.transform_to(
+ preferred_srs,
+ min(envelope.size()) / float(max(size)))
+ except core.ProjectionError:
+ # Skip projection errors
+ return None
+
+ # Calculate the zoom factor appropriate for this envelope
+ zoom = self._zoom_for_envelope(layer, pref_envelope, size)
+
+ # Get the destination raster
+ raster = _gdal.create_render_dataset(envelope, size, 4)
+ assert tuple(size) == (raster.dataset.RasterXSize, raster.dataset.RasterYSize)
+
+ # How many tiles overall in x and y at this zoom
+ n_proj_tiles = math.pow(2, zoom)
+
+ # Compute the tile co-ordinates of each corner
+ corners = [
+ layer.projection.projCoordinate(Point(pref_envelope.left, pref_envelope.top)),
+ layer.projection.projCoordinate(Point(pref_envelope.right, pref_envelope.top)),
+ layer.projection.projCoordinate(Point(pref_envelope.left, pref_envelope.bottom)),
+ layer.projection.projCoordinate(Point(pref_envelope.right, pref_envelope.bottom)),
+ ]
+ corners = [c.zoomTo(zoom) for c in corners]
+
+ corner_rows = [int(math.floor(x.row)) for x in corners]
+ corner_columns = [int(math.floor(x.column)) for x in corners]
+
+ tile_rasters = []
+
+ # Get each tile image
+ png_driver = gdal.GetDriverByName('PNG')
+ desired_srs_wkt = envelope.spatial_reference.ExportToWkt()
+ assert png_driver is not None
+ for r in xrange(min(corner_rows), max(corner_rows)+1):
+ if r < 0 or r >= n_proj_tiles:
+ continue
+ for c in xrange(min(corner_columns), max(corner_columns)+1):
+ tile_tl_point = layer.projection.coordinateProj(Coordinate(r,c,zoom))
+ tile_br_point = layer.projection.coordinateProj(Coordinate(r+1,c+1,zoom))
+
+ c = c % n_proj_tiles
+ if c < 0:
+ c += n_proj_tiles
+ tile_coord = Coordinate(r, c, zoom)
+
+ try:
+ tile_type, png_data = TileStache.getTile(layer, tile_coord, 'png')
+ if tile_type != 'image/png':
+ print('Did not get PNG data when fetching tile %s. Skipping.' % (tile_coord))
+ continue
+ except IOError as e:
+ print('Ignoring error fetching tile %s (%s).' % (tile_coord, e))
+ continue
+
+ gdal.FileFromMemBuffer('/vsimem/tmptile.png', png_data)
+ tile_raster = gdal.Open('/vsimem/tmptile.png')
+ tile_raster.SetProjection(preferred_srs_wkt)
+
+ xscale = (tile_br_point.x - tile_tl_point.x) / tile_raster.RasterXSize
+ yscale = (tile_br_point.y - tile_tl_point.y) / tile_raster.RasterXSize
+ tile_raster.SetGeoTransform((
+ tile_tl_point.x, xscale, 0.0,
+ tile_tl_point.y, 0.0, yscale,
+ ))
+
+ tile_rasters.append(Raster.from_dataset(tile_raster))
+
+ tile_raster = None
+ gdal.Unlink('/vsimem/tmptile.png')
+
+ output = Raster(
+ np.ma.masked_all((size[1], size[0], 4), dtype=np.float32),
+ envelope,
+ to_rgba=RgbaFromBands(
+ (
+ (RgbaFromBands.RED, 1.0/255.0),
+ (RgbaFromBands.GREEN, 1.0/255.0),
+ (RgbaFromBands.BLUE, 1.0/255.0),
+ (RgbaFromBands.ALPHA, 1.0/255.0),
+ ), False)
+ )
+ reproject_rasters(output, tile_rasters)
+
+ return output
+
+class TiledRasterFilter(object):
+ def __init__(self, render_cb, tile_size=None):
+ self.render_cb = render_cb
+ self.tile_size = tile_size
+
+ def __call__(self, envelope=None, size=None, **kwargs):
+ if envelope is None:
+ return None
+
+ if size is None:
+ size = map(int, envelope.size())
+
+ if self.tile_size is None:
+ raster = self.render_cb(envelope=envelope, size=size, **kwargs)
+ if raster is None or len(raster) == 0 or raster[0] is None:
+ return None
+ return raster[0]
+
+ tiles = []
+ tile_offsets = []
+ xscale, yscale = [x[0] / float(x[1]) for x in zip(envelope.offset(), size)]
+ for x in xrange(0, size[0], self.tile_size):
+ width = min(x + self.tile_size, size[0]) - x
+ for y in xrange(0, size[1], self.tile_size):
+ height = min(y + self.tile_size, size[1]) - y
+ tile_envelope = core.Envelope(
+ envelope.left + xscale * x,
+ envelope.left + xscale * (x + width),
+ envelope.top + yscale * y,
+ envelope.top + yscale * (y + height),
+ envelope.spatial_reference,
+ )
+ tiles.append(self.render_cb(envelope=tile_envelope, size=(width, height), **kwargs))
+ tile_offsets.append((x,y))
+
+ results = [x for x in zip(tile_offsets, tiles) if x[1] is not None]
+ if len(results) == 0:
+ return None
+
+ # FIXME: This assumes that to_rgba_cb is the same for each tile
+ prototype = results[0][1]
+ depth = prototype.array.shape[2]
+
+ shape = (size[1], size[0])
+ mask = np.ones(shape + (depth,), dtype=np.bool)
+ data = np.zeros(shape + (depth,), dtype=np.float32)
+ for pos, tile in results:
+ x, y = pos
+ h, w = tile.array.shape[:2]
+
+ tile_mask = np.ma.getmask(tile.array)
+ if tile_mask is np.ma.nomask:
+ tile_mask = False
+
+ mask[y:(y+h), x:(x+w), :] = tile_mask
+ data[y:(y+h), x:(x+w), :] = tile.array
+
+ if np.any(mask):
+ data = np.ma.array(data, mask=mask)
+
+ return Raster(data, envelope, prototype=prototype)
+
+class ReprojectingRasterFilter(object):
+ def __init__(self, native_spatial_reference, render_cb):
+ self.native_spatial_reference = native_spatial_reference
+ self.native_spatial_reference_wkt = native_spatial_reference.ExportToWkt()
+ self.render_cb = render_cb
+
+ def __call__(self, envelope=None, size=None, **kwargs):
+ if envelope is None:
+ return None
+
+ if size is None:
+ size = map(int, envelope.size())
+
+ if envelope.spatial_reference.IsSame(self.native_spatial_reference):
+ return self.render_cb(envelope=envelope, size=size, **kwargs)
+
+ # We need to reproject this data. Convert the envelope into the native spatial reference
+ try:
+ native_envelope = envelope.transform_to(
+ self.native_spatial_reference,
+ min(envelope.size()) / float(max(size)))
+ except ProjectionError:
+ # If we fail to reproject, return a null tile
+ return None
+
+ # Get the native tile
+ raster = self.render_cb(envelope=native_envelope, size=size, **kwargs)
+ if raster is None:
+ return None
+
+ raster_ds = raster.as_dataset()
+
+ # Get the destination raster
+ ds = _gdal.create_render_dataset(envelope, size, 4).dataset
+
+ desired_srs_wkt = envelope.spatial_reference.ExportToWkt()
+ gdal.ReprojectImage(
+ raster_ds, ds,
+ self.native_spatial_reference_wkt,
+ desired_srs_wkt,
+ gdal.GRA_Bilinear if raster.can_interpolate else gdal.GRA_NearestNeighbour)
+
+ # Create a mask raster
+ mask_raster = _gdal.create_render_dataset(envelope, size, 1, data_type=gdal.GDT_Float32)
+ mask_ds = mask_raster.dataset
+ band = mask_ds.GetRasterBand(1)
+ band.SetColorInterpretation(gdal.GCI_Undefined)
+ band.SetNoDataValue(float('nan'))
+ band.Fill(float('nan'))
+ gdal.ReprojectImage(
+ raster_ds, mask_ds,
+ self.native_spatial_reference_wkt,
+ desired_srs_wkt,
+ gdal.GRA_NearestNeighbour)
+ mask_band = mask_ds.GetRasterBand(1).GetMaskBand()
+
+ return Raster.from_dataset(ds, mask_band=mask_band, prototype=raster)
View
163 src/foldbeam/tests.py
@@ -1,12 +1,14 @@
-from . import core, graph, nodes, transform, pads, _gdal
-from . graph import connect
import math
+import os
+import unittest
+
import numpy as np
from osgeo import gdal
from osgeo.osr import SpatialReference
from TileStache.Config import buildConfiguration
-import os
-import unittest
+
+from foldbeam import core, graph, _gdal
+from foldbeam import raster, vector, tilestache
class TestUtility(unittest.TestCase):
def test_create_render_dataset(self):
@@ -50,95 +52,94 @@ def test_tilestache_lat_lng(self):
envelope_srs.ImportFromEPSG(4326) # WGS 84 lat/lng
envelope = core.Envelope(-180, 180, 89, -89, envelope_srs)
- node = nodes.TileStacheRasterNode(self.config.layers['osm'])
+ node = raster.TileStacheSource(config=self.config)
size = (1024, 512)
- raster = node.outputs.output(envelope=envelope, size=size)
- self.assertIsNotNone(raster)
- self.assertIsInstance(raster, core.Raster)
- self.assertEqual(raster.array.shape[1], size[0])
- self.assertEqual(raster.array.shape[0], size[1])
+ tile = node.outputs.osm(envelope=envelope, size=size)
+ self.assertIsNotNone(tile)
+ self.assertIsInstance(tile, raster.Raster)
+ self.assertEqual(tile.array.shape[1], size[0])
+ self.assertEqual(tile.array.shape[0], size[1])
- ds = raster.as_dataset()
+ ds = tile.as_dataset()
self.assertEqual(ds.RasterXSize, size[0])
self.assertEqual(ds.RasterYSize, size[1])
- raster.write_tiff('world-latlng.tiff')
+ tile.write_tiff('world-latlng.tiff')
def test_tilestache_usna(self):
envelope_srs = SpatialReference()
envelope_srs.ImportFromEPSG(2163) # US National Atlas Equal Area
envelope = core.Envelope(-3e6, 3e6, 2e6, -2e6, envelope_srs)
- node = nodes.TileStacheRasterNode(self.config.layers['osm'])
+ node = raster.TileStacheSource(config=self.config)
size = (1200, 800)
- raster = node.outputs.output(envelope=envelope, size=size)
- self.assertIsNotNone(raster)
- self.assertIsInstance(raster, core.Raster)
- self.assertEqual(raster.array.shape[1], size[0])
- self.assertEqual(raster.array.shape[0], size[1])
+ tile = node.outputs.osm(envelope=envelope, size=size)
+ self.assertIsNotNone(tile)
+ self.assertIsInstance(tile, raster.Raster)
+ self.assertEqual(tile.array.shape[1], size[0])
+ self.assertEqual(tile.array.shape[0], size[1])
- ds = raster.as_dataset()
+ ds = tile.as_dataset()
self.assertEqual(ds.RasterXSize, size[0])
self.assertEqual(ds.RasterYSize, size[1])
- raster.write_tiff('world-usna.tiff')
+ tile.write_tiff('world-usna.tiff')
def test_tilestache_osgrid_overlay(self):
envelope_srs = SpatialReference()
envelope_srs.ImportFromEPSG(27700) # British national grid
envelope = core.Envelope(0, 700000, 1300000, 0, envelope_srs)
- roads_node = nodes.TileStacheRasterNode(self.config.layers['osm'])
- aerial_node = nodes.TileStacheRasterNode(self.config.layers['aerial'])
- node = nodes.LayerRasterNode(top_opacity=0.5)
- connect(roads_node.outputs.output, node.inputs.top)
- connect(aerial_node.outputs.output, node.inputs.bottom)
+ ts_node = raster.TileStacheSource(config=self.config)
+ node = raster.CompositeOver(top_opacity=0.5)
+ graph.connect(ts_node.outputs.osm, node.inputs.top)
+ graph.connect(ts_node.outputs.aerial, node.inputs.bottom)
size = (700, 1300)
- raster = node.outputs.output(envelope=envelope, size=size)
- self.assertIsNotNone(raster)
- self.assertIsInstance(raster, core.Raster)
- self.assertEqual(raster.array.shape[1], size[0])
- self.assertEqual(raster.array.shape[0], size[1])
+ tile = node.outputs.output(envelope=envelope, size=size)
+ self.assertIsNotNone(tile)
+ self.assertIsInstance(tile, raster.Raster)
+ self.assertEqual(tile.array.shape[1], size[0])
+ self.assertEqual(tile.array.shape[0], size[1])
- ds = raster.as_dataset()
+ ds = tile.as_dataset()
self.assertEqual(ds.RasterXSize, size[0])
self.assertEqual(ds.RasterYSize, size[1])
- raster.write_tiff('uk-osgrid-overlay.tiff')
+ tile.write_tiff('uk-osgrid-overlay.tiff')
def test_tilestache_osgrid(self):
envelope_srs = SpatialReference()
envelope_srs.ImportFromEPSG(27700) # British national grid
envelope = core.Envelope(0, 700000, 1300000, 0, envelope_srs)
- node = nodes.TileStacheRasterNode(self.config.layers['osm'])
+ node = raster.TileStacheSource(config=self.config)
size = (700, 1300)
- raster = node.outputs.output(envelope=envelope, size=size)
- self.assertIsNotNone(raster)
- self.assertIsInstance(raster, core.Raster)
- self.assertEqual(raster.array.shape[1], size[0])
- self.assertEqual(raster.array.shape[0], size[1])
+ tile = node.outputs.osm(envelope=envelope, size=size)
+ self.assertIsNotNone(tile)
+ self.assertIsInstance(tile, raster.Raster)
+ self.assertEqual(tile.array.shape[1], size[0])
+ self.assertEqual(tile.array.shape[0], size[1])
- ds = raster.as_dataset()
+ ds = tile.as_dataset()
self.assertEqual(ds.RasterXSize, size[0])
self.assertEqual(ds.RasterYSize, size[1])
- raster.write_tiff('uk-osgrid.tiff')
+ tile.write_tiff('uk-osgrid.tiff')
def test_tilestache_osgrid_crazy(self):
envelope_srs = SpatialReference()
envelope_srs.ImportFromEPSG(27700) # British national grid
envelope = core.Envelope(-2e6, 2e6, 3e6, -3e6, envelope_srs)
- node = nodes.TileStacheRasterNode(self.config.layers['osm'])
+ node = raster.TileStacheSource(config=self.config)
size = (800, 1200)
- raster = node.outputs.output(envelope=envelope, size=size)
- self.assertIsNotNone(raster)
- self.assertIsInstance(raster, core.Raster)
- self.assertEqual(raster.array.shape[1], size[0])
- self.assertEqual(raster.array.shape[0], size[1])
+ tile = node.outputs.osm(envelope=envelope, size=size)
+ self.assertIsNotNone(tile)
+ self.assertIsInstance(tile, raster.Raster)
+ self.assertEqual(tile.array.shape[1], size[0])
+ self.assertEqual(tile.array.shape[0], size[1])
- ds = raster.as_dataset()
+ ds = tile.as_dataset()
self.assertEqual(ds.RasterXSize, size[0])
self.assertEqual(ds.RasterYSize, size[1])
- raster.write_tiff('mad-osgrid.tiff')
+ tile.write_tiff('mad-osgrid.tiff')
def test_tilestache_big_ben(self):
# square around Big Ben
@@ -151,18 +152,18 @@ def test_tilestache_big_ben(self):
centre[0]-0.5*skirt, centre[0]+0.5*skirt, centre[1]+0.5*skirt, centre[1]-0.5*skirt,
envelope_srs)
- node = nodes.TileStacheRasterNode(self.config.layers['osm'])
+ node = raster.TileStacheSource(config=self.config)
size = (512, 512)
- raster = node.outputs.output(envelope=envelope, size=size)
- self.assertIsNotNone(raster)
- self.assertIsInstance(raster, core.Raster)
- self.assertEqual(raster.array.shape[1], size[0])
- self.assertEqual(raster.array.shape[0], size[1])
+ tile = node.outputs.osm(envelope=envelope, size=size)
+ self.assertIsNotNone(tile)
+ self.assertIsInstance(tile, raster.Raster)
+ self.assertEqual(tile.array.shape[1], size[0])
+ self.assertEqual(tile.array.shape[0], size[1])
- ds = raster.as_dataset()
+ ds = tile.as_dataset()
self.assertEqual(ds.RasterXSize, size[0])
self.assertEqual(ds.RasterYSize, size[1])
- raster.write_tiff('big-ben.tiff')
+ tile.write_tiff('big-ben.tiff')
def test_tilestache_big_ben_small(self):
# square around Big Ben
@@ -175,18 +176,18 @@ def test_tilestache_big_ben_small(self):
centre[0]-0.5*skirt, centre[0]+0.5*skirt, centre[1]+0.5*skirt, centre[1]-0.5*skirt,
envelope_srs)
- node = nodes.TileStacheRasterNode(self.config.layers['osm'])
+ node = raster.TileStacheSource(config=self.config)
size = (512, 512)
- raster = node.outputs.output(envelope=envelope, size=size)
- self.assertIsNotNone(raster)
- self.assertIsInstance(raster, core.Raster)
- self.assertEqual(raster.array.shape[1], size[0])
- self.assertEqual(raster.array.shape[0], size[1])
+ tile = node.outputs.osm(envelope=envelope, size=size)
+ self.assertIsNotNone(tile)
+ self.assertIsInstance(tile, raster.Raster)
+ self.assertEqual(tile.array.shape[1], size[0])
+ self.assertEqual(tile.array.shape[0], size[1])
- ds = raster.as_dataset()
+ ds = tile.as_dataset()
self.assertEqual(ds.RasterXSize, size[0])
self.assertEqual(ds.RasterYSize, size[1])
- raster.write_tiff('big-ben-sm.tiff')
+ tile.write_tiff('big-ben-sm.tiff')
def test_tilestache_big_ben_proj_units(self):
# test letting image size be determined from projection units
@@ -201,17 +202,17 @@ def test_tilestache_big_ben_proj_units(self):
centre[0]-0.5*skirt[0], centre[0]+0.5*skirt[0], centre[1]+0.5*skirt[1], centre[1]-0.5*skirt[1],
envelope_srs)
- node = nodes.TileStacheRasterNode(self.config.layers['osm'])
- raster = node.outputs.output(envelope=envelope)
- self.assertIsNotNone(raster)
- self.assertIsInstance(raster, core.Raster)
- self.assertEqual(raster.array.shape[1], skirt[0])
- self.assertEqual(raster.array.shape[0], skirt[1])
+ node = raster.TileStacheSource(config=self.config)
+ tile = node.outputs.osm(envelope=envelope)
+ self.assertIsNotNone(tile)
+ self.assertIsInstance(tile, raster.Raster)
+ self.assertEqual(tile.array.shape[1], skirt[0])
+ self.assertEqual(tile.array.shape[0], skirt[1])
- ds = raster.as_dataset()
+ ds = tile.as_dataset()
self.assertEqual(ds.RasterXSize, skirt[0])
self.assertEqual(ds.RasterYSize, skirt[1])
- raster.write_tiff('big-ben-metres.tiff')
+ tile.write_tiff('big-ben-metres.tiff')
class TestBoundary(unittest.TestCase):
def test_bbox(self):
@@ -242,7 +243,7 @@ def test_damage(self):
def damage_cb(d):
damages.append(d)
- e = graph.OutputPad(graph.RasterType, self, 'foo', lambda: None)
+ e = graph.OutputPad(raster.Raster, self, 'foo', lambda: None)
e.damaged.connect(damage_cb)
self.assertEqual(len(damages), 0)
@@ -282,7 +283,7 @@ def setUp(self):
lnglat_srs)
def test_compute_warp(self):
- ll_x, ll_y = transform.compute_warp(
+ ll_x, ll_y = raster.compute_warp(
self.bng_env, (256, 256),
self.lnglat_env.spatial_reference)
self.assertTrue(np.all(ll_x > -2))
@@ -293,7 +294,7 @@ def test_compute_warp(self):
def test_proj_to_pixels(self):
x = np.array([[-1.4, 0.1], [-2, 2], [-3, 8]])
y = np.array([[50.2, 50.4], [54, 50], [59, 49]])
- px, py = transform.proj_to_pixels(x, y, self.lnglat_env, (128, 64))
+ px, py = raster.proj_to_pixels(x, y, self.lnglat_env, (128, 64))
self.assertGreaterEqual(px[0,0], 0)
self.assertLessEqual(px[0,0], 128)
@@ -318,7 +319,7 @@ def test_proj_to_pixels(self):
self.assertGreater(py[2,0], 64)
self.assertLess(py[2,1], 0)
- ll_x, ll_y = transform.compute_warp(
+ ll_x, ll_y = raster.compute_warp(
self.bng_env, (256, 256),
self.lnglat_env.spatial_reference)
self.assertTrue(np.all(ll_x > -2))
@@ -329,7 +330,7 @@ def test_proj_to_pixels(self):
def test_sample_raster(self):
src_data = np.atleast_3d(np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]]))
src_data = np.dstack((src_data, 2*src_data))
- src = core.Raster(src_data, self.bng_env)
+ src = raster.Raster(src_data, self.bng_env)
src_x = np.array([
[-1.2, -0.2, 0.8, 1.8, 2.8, 3.2, 4.2],
[-1.3, -0.3, 0.7, 1.7, 2.7, 3.2, 4.2],
@@ -342,7 +343,7 @@ def test_sample_raster(self):
[0, 1, 2, 2, 2.3, 1.2, 0.2],
[0, 1, 2, 2, 2.3, 1.2, 0.2],
])
- dst = transform.sample_raster(src_x, src_y, src)
+ dst = raster.sample_raster(src_x, src_y, src)
self.assertEqual(dst.shape[:2], src_x.shape[:2])
self.assertEqual(dst.shape[2], src_data.shape[2])
@@ -356,10 +357,10 @@ def test_sample_raster(self):
def test_reproject_raster(self):
src_data = np.atleast_3d(np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]]))
src_data = np.dstack((src_data, 2*src_data))
- src = core.Raster(src_data, self.bng_env)
+ src = raster.Raster(src_data, self.bng_env)
- dst = core.Raster(np.ones((10,5)), self.lnglat_env)
- transform.reproject_raster(dst, src)
+ dst = raster.Raster(np.ones((10,5)), self.lnglat_env)
+ raster.reproject_raster(dst, src)
def test_suite():
return unittest.TestSuite([
View
21 src/foldbeam/tilestache.py
@@ -1,18 +1,23 @@
-from .core import Envelope
-from .pipeline import Pipeline
-from .graph import Node, node, RasterType
+"""
+Serving output via TileStache
+=============================
+
+"""
+
import json
+
import numpy as np
from osgeo import osr
from PIL import Image
import TileStache
-from werkzeug.serving import run_simple
-@node
-class TileStacheServerNode(Node):
+from foldbeam import core, graph, raster
+
+@graph.node
+class TileStacheServerNode(graph.Node):
def __init__(self):
super(TileStacheServerNode, self).__init__()
- self.add_input('raster', RasterType)
+ self.add_input('raster', raster.Raster)
# Create the tilestache config
self.config = TileStache.Config.buildConfiguration({
@@ -38,7 +43,7 @@ def __init__(self, raster_pad):
def renderArea(self, width, height, srs, xmin, ymin, xmax, ymax, zoom):
spatial_reference = osr.SpatialReference()
spatial_reference.ImportFromProj4(srs)
- envelope = Envelope(xmin, xmax, ymax, ymin, spatial_reference)
+ envelope = core.Envelope(xmin, xmax, ymax, ymin, spatial_reference)
raster = self.raster_pad(envelope=envelope, size=(width, height))
if raster is None:
View
9 src/foldbeam/tool/render.py
@@ -1,6 +1,6 @@
import argparse
import TileStache
-from foldbeam import nodes, core, graph, pads
+from foldbeam import core, graph, raster
from osgeo import gdal
from osgeo.osr import SpatialReference
import sys
@@ -81,10 +81,11 @@ def run(args):
ew, eh = map(abs, envelope.offset())
args.width = max(1, int(args.height * ew / eh))
- node = nodes.TileStacheRasterNode(config.layers['aerial' if args.aerial else 'osm'])
+ node = raster.TileStacheSource(config=config)
size = (args.width, args.height)
- raster = node.outputs.output(envelope=envelope, size=size)
- raster.write_tiff(args.output)
+ output_pad = node.outputs['aerial' if args.aerial else 'osm']
+ output = output_pad(envelope=envelope, size=size)
+ output.write_tiff(args.output)
if __name__ == '__main__':
main()
View
171 src/foldbeam/transform.py
@@ -1,171 +0,0 @@
-"""
-Geographic co-ordinate transformation
-=====================================
-
-This module contains support functions for transforming between different geographic co-ordinate systems.
-
-"""
-
-from . import core
-import numpy as np
-import os
-from osgeo import gdal
-import pyproj
-
-def compute_warp(dst_envelope, dst_size, src_spatial_reference):
- """Compute the projection co-ordinates in a source spatial projection for pixels in a destination image.
-
- :param dst_envelope: the envelope of the destination raster
- :type dst_envelope: :py:class:`core.Envelope`
- :param dst_size: the width and height of the destination raster
- :type dst_size: pair of integer
- :param src_spatial_reference: co-ordinate system of the source raster
- :type src_spatial_reference: :py:class:`osgeo.osr.SpatialReference`
- :rtype: pair of arrays
-
- Compute the projection co-ordinates of each pixel in a destination raster when projected into the source spatial
- reference. This is necessary to re-project images correctly.
-
- The return value is a pair of array-like objects giving the x and y co-ordinates. The returned arrays are the width
- and height specified in dst_size.
-
- """
-
- # Compute x- and y-co-ordinates in dst from the envelope
- dst_cols, dst_rows = dst_size
- dst_x = np.repeat(
- np.linspace(dst_envelope.left, dst_envelope.right, dst_cols).reshape(1, dst_cols),
- dst_rows, 0)
- dst_y = np.repeat(
- np.linspace(dst_envelope.top, dst_envelope.bottom, dst_rows).reshape(dst_rows, 1),
- dst_cols, 1)
-
- assert dst_x.shape == (dst_rows, dst_cols)
- assert dst_y.shape == (dst_rows, dst_cols)
-
- # Special case when source and destination have same srs
- if dst_envelope.spatial_reference.IsSame(src_spatial_reference):
- return dst_x, dst_y
-
- dst_proj = pyproj.Proj(dst_envelope.spatial_reference.ExportToProj4())
- src_proj = pyproj.Proj(src_spatial_reference.ExportToProj4())
-
- src_x, src_y = pyproj.transform(dst_proj, src_proj, dst_x, dst_y)
-
- return src_x, src_y
-
-def proj_to_pixels(x, y, envelope, size):
- """Compute pixel co-ordinates of projection co-ordinates x and y in image with envelope and size specified."""
-
- off = envelope.offset()
- x = (x - envelope.left) * (float(size[0]-1) / off[0])
- y = (y - envelope.top) * (float(size[1]-1) / off[1])
- return x,y
-
-def pixels_to_proj(x, y, envelope, size):
- """Compute projection co-ordinates of pixel co-ordinates x and y in image with envelope and size specified."""
-
- off = envelope.offset()
- x = x * (off[0] / float(size[0]-1)) + envelope.left
- y = y * (off[1] / float(size[1]-1)) + envelope.right
- return x,y
-
-def sample_raster(src_pixel_x, src_pixel_y, src_raster):
- """Sample pixels from a raster.
-
- :param src_pixel_x: array of pixel x-co-ordinates
- :param src_pixel_y: array of pixel y-co-ordinates
- :param src_raster: source raster
- :type src_raster: :py:class:`core.Raster`
- :rtype: :py:class:`numpy.array` or :py:class:`numpy.ma.array`
-
- The returned array has the same shape as src_raster with each value corresponding to the sampled value. The arrays
- src_pixel_x and src_pixel_y must have the same width and height as src_raster. If pixel co-ordinates outside of the
- range of those present in src_raster are specified, the values will be masked out of the resulting array and a numpy
- masked array will be returned.
-
- """
-
- px = np.int32(np.round(src_pixel_x))
- py = np.int32(np.round(src_pixel_y))
-
- x_valid = np.logical_and(px >= 0, px < src_raster.array.shape[1])
- y_valid = np.logical_and(py >= 0, py < src_raster.array.shape[0])
- valid = np.logical_and(x_valid, y_valid)
-
- valid_flag = valid.ravel()
- valid_src_indices = np.ravel_multi_index(
- (py.ravel()[valid_flag], px.ravel()[valid_flag]),
- src_raster.array.shape[:2]
- )
-
- dst_shape = np.atleast_2d(src_pixel_x).shape[:2] + (src_raster.array.shape[2],)
- dst = np.empty(shape=dst_shape, dtype=np.float32)
-
- for i in xrange(dst_shape[2]):
- dst[valid, i] = src_raster.array[:,:,i].ravel()[valid_src_indices]
-
- dst = np.ma.fix_invalid(
- dst,
- fill_value=np.nan,
- mask=np.repeat(np.atleast_3d(np.logical_not(valid)), dst_shape[2], 2)
- )
-
- return dst
-
-def reproject_raster(dst, src):
- """Re-project a source raster into a destination raster.
-
- See :py:func:`reproject_rasters`.
-
- """
- reproject_rasters(dst, (src,))
-
-def _py_reproject_rasters(dst, srcs):
- if len(srcs) == 0:
- return
-
- src_srs = srcs[0].envelope.spatial_reference
- assert all([x.en