Newer
Older
# -*- coding: utf-8 -*-
"""
/***************************************************************************
HUB TimeSeriesViewer
-------------------
begin : 2015-08-20
git sha : $Format:%H$
copyright : (C) 2017 by HU-Berlin
email : benjamin.jakimow@geo.hu-berlin.de
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
"""
# noinspection PyPep8Naming
from __future__ import absolute_import

benjamin.jakimow@geo.hu-berlin.de
committed
import os, sys, re, pickle, tempfile
from collections import OrderedDict
import tempfile

benjamin.jakimow@geo.hu-berlin.de
committed
from osgeo import gdal, osr, ogr

benjamin.jakimow@geo.hu-berlin.de
committed
from qgis.gui import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *

benjamin.jakimow@geo.hu-berlin.de
committed

benjamin.jakimow@geo.hu-berlin.de
committed
def px2geo(px, gt):
#see http://www.gdal.org/gdal_datamodel.html
gx = gt[0] + px.x()*gt[1]+px.y()*gt[2]
gy = gt[3] + px.x()*gt[4]+px.y()*gt[5]
return QgsPoint(gx,gy)
class VRTRasterInputSourceBand(object):
def __init__(self, path, bandIndex, bandName=''):
self.mPath = os.path.normpath(path)
self.mBandIndex = bandIndex
self.mBandName = bandName
self.mNoData = None

benjamin.jakimow@geo.hu-berlin.de
committed

benjamin.jakimow@geo.hu-berlin.de
committed
def isEqual(self, other):
if isinstance(other, VRTRasterInputSourceBand):
return self.mPath == other.mPath and self.mBandIndex == other.mBandIndex

benjamin.jakimow@geo.hu-berlin.de
committed
def __reduce_ex__(self, protocol):
return self.__class__, (self.mPath, self.mBandIndex, self.mBandName), self.__getstate__()
def __getstate__(self):
state = self.__dict__.copy()
state.pop('mVirtualBand')
return state
def __setstate__(self, state):
self.__dict__.update(state)
def virtualBand(self):
return self.mVirtualBand

benjamin.jakimow@geo.hu-berlin.de
committed
class VRTRasterBand(QObject):
sigNameChanged = pyqtSignal(str)
sigSourceInserted = pyqtSignal(int, VRTRasterInputSourceBand)
sigSourceRemoved = pyqtSignal(int, VRTRasterInputSourceBand)
def __init__(self, name='', parent=None):

benjamin.jakimow@geo.hu-berlin.de
committed
super(VRTRasterBand, self).__init__(parent)
self.sources = []
self.mName = name

benjamin.jakimow@geo.hu-berlin.de
committed
def setName(self, name):
oldName = self.mName
self.mName = name
if oldName != self.mName:
self.sigNameChanged.emit(name)
def name(self):
return self.mName

benjamin.jakimow@geo.hu-berlin.de
committed
def addSource(self, virtualBandInputSource):
assert isinstance(virtualBandInputSource, VRTRasterInputSourceBand)
self.insertSource(len(self.sources), virtualBandInputSource)
def insertSource(self, index, virtualBandInputSource):
assert isinstance(virtualBandInputSource, VRTRasterInputSourceBand)
virtualBandInputSource.mVirtualBand = self
assert index <= len(self.sources)
self.sources.insert(index, virtualBandInputSource)

benjamin.jakimow@geo.hu-berlin.de
committed
self.sigSourceInserted.emit(index, virtualBandInputSource)

benjamin.jakimow@geo.hu-berlin.de
committed
if isinstance(self.mVRT, VRTRaster):
return self.mVRT.mBands.index(self)

benjamin.jakimow@geo.hu-berlin.de
committed
def removeSource(self, vrtRasterInputSourceBand):

benjamin.jakimow@geo.hu-berlin.de
committed
Removes a VRTRasterInputSourceBand
:param vrtRasterInputSourceBand: band index| VRTRasterInputSourceBand
:return: The VRTRasterInputSourceBand that was removed

benjamin.jakimow@geo.hu-berlin.de
committed
if not isinstance(vrtRasterInputSourceBand, VRTRasterInputSourceBand):
vrtRasterInputSourceBand = self.sources[vrtRasterInputSourceBand]
if vrtRasterInputSourceBand in self.sources:
i = self.sources.index(vrtRasterInputSourceBand)
self.sources.remove(vrtRasterInputSourceBand)
self.sigSourceRemoved.emit(i, vrtRasterInputSourceBand)
def sourceFiles(self):
"""
:return: list of file-paths to all source files
"""

benjamin.jakimow@geo.hu-berlin.de
committed
files = set([inputSource.mPath for inputSource in self.sources])
return sorted(list(files))
def __repr__(self):
infos = ['VirtualBand name="{}"'.format(self.mName)]
for i, info in enumerate(self.sources):

benjamin.jakimow@geo.hu-berlin.de
committed
assert isinstance(info, VRTRasterInputSourceBand)
infos.append('\t{} SourceFileName {} SourceBand {}'.format(i + 1, info.mPath, info.mBandIndex))
return '\n'.join(infos)

benjamin.jakimow@geo.hu-berlin.de
committed
LUT_ReampleAlg = {'nearest': gdal.GRA_NearestNeighbour,
'bilinear': gdal.GRA_Bilinear,
'mode':gdal.GRA_Mode,
'lanczos':gdal.GRA_Lanczos,
'average':gdal.GRA_Average,
'cubic':gdal.GRA_Cubic,
'cubic_splie':gdal.GRA_CubicSpline}
class VRTRasterPreviewMapCanvas(QgsMapCanvas):
def __init__(self, parent=None, *args, **kwds):
super(VRTRasterPreviewMapCanvas, self).__init__(parent, *args, **kwds)

benjamin.jakimow@geo.hu-berlin.de
committed
self.setCrsTransformEnabled(True)
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
def contextMenuEvent(self, event):
menu = QMenu()
action = menu.addAction('Refresh')
action.triggered.connect(self.refresh)
action = menu.addAction('Reset')
action.triggered.connect(self.reset)
menu.exec_(event.globalPos())
def setLayerSet(self, layers):
raise DeprecationWarning()
def setLayers(self, layers):
assert isinstance(layers, list)
def area(layer):
extent = layer.extent()
return extent.width() * extent.height()
layers = list(sorted(layers, key = lambda lyr: area(lyr), reverse=True))
QgsMapLayerRegistry.instance().addMapLayers(layers)
super(VRTRasterPreviewMapCanvas, self).setLayerSet([QgsMapCanvasLayer(l) for l in layers])
def reset(self):
extent = self.fullExtent()
extent.scale(1.05)
self.setExtent(extent)
self.refresh()

benjamin.jakimow@geo.hu-berlin.de
committed
class VRTRaster(QObject):

benjamin.jakimow@geo.hu-berlin.de
committed
sigSourceBandInserted = pyqtSignal(VRTRasterBand, VRTRasterInputSourceBand)
sigSourceBandRemoved = pyqtSignal(VRTRasterBand, VRTRasterInputSourceBand)
sigSourceRasterAdded = pyqtSignal(list)
sigSourceRasterRemoved = pyqtSignal(list)
sigBandInserted = pyqtSignal(int, VRTRasterBand)
sigBandRemoved = pyqtSignal(int, VRTRasterBand)
sigCrsChanged = pyqtSignal(QgsCoordinateReferenceSystem)

benjamin.jakimow@geo.hu-berlin.de
committed
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
def __init__(self, parent=None):
super(VRTRaster, self).__init__(parent)
self.mBands = []
self.mCrs = None
self.mResampleAlg = gdal.GRA_NearestNeighbour
self.mMetadata = dict()
self.mSourceRasterBounds = dict()
self.mOutputBounds = None
self.sigSourceBandRemoved.connect(self.updateSourceRasterBounds)
self.sigSourceBandInserted.connect(self.updateSourceRasterBounds)
self.sigBandRemoved.connect(self.updateSourceRasterBounds)
self.sigBandInserted.connect(self.updateSourceRasterBounds)
def setCrs(self, crs):
if isinstance(crs, osr.SpatialReference):
auth = '{}:{}'.format(crs.GetAttrValue('AUTHORITY',0), crs.GetAttrValue('AUTHORITY',1))
crs = QgsCoordinateReferenceSystem(auth)
if isinstance(crs, QgsCoordinateReferenceSystem):
if crs != self.mCrs:
self.mCrs = crs
self.sigCrsChanged.emit(self.mCrs)
def crs(self):
return self.mCrs
def addVirtualBand(self, virtualBand):
"""
Adds a virtual band
:param virtualBand: the VirtualBand to be added
:return: VirtualBand
"""

benjamin.jakimow@geo.hu-berlin.de
committed
assert isinstance(virtualBand, VRTRasterBand)
return self.insertVirtualBand(len(self), virtualBand)
def insertSourceBand(self, virtualBandIndex, pathSource, sourceBandIndex):
"""
Inserts a source band into the VRT stack
:param virtualBandIndex: target virtual band index
:param pathSource: path of source file
:param sourceBandIndex: source file band index
"""

benjamin.jakimow@geo.hu-berlin.de
committed
while virtualBandIndex > len(self.mBands)-1:
self.insertVirtualBand(len(self.mBands), VRTRasterBand())

benjamin.jakimow@geo.hu-berlin.de
committed
vBand = self.mBands[virtualBandIndex]
vBand.addSourceBand(pathSource, sourceBandIndex)

benjamin.jakimow@geo.hu-berlin.de
committed
def insertVirtualBand(self, index, virtualBand):

benjamin.jakimow@geo.hu-berlin.de
committed
:param index: the insert position
:param virtualBand: the VirtualBand to be inserted
:return: the VirtualBand
"""

benjamin.jakimow@geo.hu-berlin.de
committed
assert isinstance(virtualBand, VRTRasterBand)
assert index <= len(self.mBands)
if len(virtualBand.name()) == 0:
virtualBand.setName('Band {}'.format(index+1))

benjamin.jakimow@geo.hu-berlin.de
committed
virtualBand.sigSourceInserted.connect(
lambda _, sourceBand: self.sigSourceBandInserted.emit(virtualBand, sourceBand))
virtualBand.sigSourceRemoved.connect(
lambda _, sourceBand: self.sigSourceBandInserted.emit(virtualBand, sourceBand))
self.mBands.insert(index, virtualBand)
self.sigBandInserted.emit(index, virtualBand)
return self[index]
def removeVirtualBands(self, bandsOrIndices):
assert isinstance(bandsOrIndices, list)
to_remove = []

benjamin.jakimow@geo.hu-berlin.de
committed
for virtualBand in bandsOrIndices:
if not isinstance(virtualBand, VRTRasterBand):
virtualBand = self.mBands[virtualBand]
to_remove.append((self.mBands.index(virtualBand), virtualBand))

benjamin.jakimow@geo.hu-berlin.de
committed
to_remove = sorted(to_remove, key=lambda t: t[0], reverse=True)
for index, virtualBand in to_remove:
self.mBands.remove(virtualBand)
self.sigBandRemoved.emit(index, virtualBand)
def removeInputSource(self, path):
assert path in self.sourceRaster()
for vBand in self.mBands:
assert isinstance(vBand, VRTRasterBand)
if path in vBand.sources():
vBand.removeSource(path)

benjamin.jakimow@geo.hu-berlin.de
committed
self.removeVirtualBands([bandOrIndex])
def addFilesAsMosaic(self, files):
"""
Shortcut to mosaic all input files. All bands will maintain their band position in the virtual file.
:param files: [list-of-file-paths]
"""
for file in files:
ds = gdal.Open(file)
assert isinstance(ds, gdal.Dataset)
nb = ds.RasterCount
for b in range(nb):
if b+1 < len(self):
#add new virtual band

benjamin.jakimow@geo.hu-berlin.de
committed
self.addVirtualBand(VRTRasterBand())

benjamin.jakimow@geo.hu-berlin.de
committed
assert isinstance(vBand, VRTRasterBand)
def addFilesAsStack(self, files):
"""
Shortcut to stack all input files, i.e. each band of an input file will be a new virtual band.
Bands in the virtual file will be ordered as file1-band1, file1-band n, file2-band1, file2-band,...
:param files: [list-of-file-paths]
:return: self
"""
for file in files:
ds = gdal.Open(file)
assert isinstance(ds, gdal.Dataset), 'Can not open {}'.format(file)
nb = ds.RasterCount
ds = None
for b in range(nb):
#each new band is a new virtual band

benjamin.jakimow@geo.hu-berlin.de
committed
vBand = self.addVirtualBand(VRTRasterBand())
assert isinstance(vBand, VRTRasterBand)
vBand.addSource(VRTRasterInputSourceBand(file, b))

benjamin.jakimow@geo.hu-berlin.de
committed
def sourceRaster(self):
files = set()

benjamin.jakimow@geo.hu-berlin.de
committed
for vBand in self.mBands:
assert isinstance(vBand, VRTRasterBand)
files.update(set(vBand.sourceFiles()))
return sorted(list(files))

benjamin.jakimow@geo.hu-berlin.de
committed
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
def sourceRasterBounds(self):
return self.mSourceRasterBounds
def outputBounds(self):
if isinstance(self.mOutputBounds, RasterBounds):
return
#calculate from source rasters
def setOutputBounds(self, bounds):
assert isinstance(self, RasterBounds)
self.mOutputBounds = bounds
def updateSourceRasterBounds(self):
srcFiles = self.sourceRaster()
toRemove = [f for f in self.mSourceRasterBounds.keys() if f not in srcFiles]
toAdd = [f for f in srcFiles if f not in self.mSourceRasterBounds.keys()]
for f in toRemove:
del self.mSourceRasterBounds[f]
for f in toAdd:
self.mSourceRasterBounds[f] = RasterBounds(f)
if len(srcFiles) > 0 and self.crs() == None:
self.setCrs(self.mSourceRasterBounds[srcFiles[0]].crs)
elif len(srcFiles) == 0:
self.setCrs(None)
if len(toRemove) > 0:
self.sigSourceRasterRemoved.emit(toRemove)
if len(toAdd) > 0:
self.sigSourceRasterAdded.emit(toAdd)
def saveVRT(self, pathVRT, resampleAlg=gdal.GRA_NearestNeighbour, **kwds):
"""
:param pathVRT: path to VRT that is created
:param options --- can be be an array of strings, a string or let empty and filled from other keywords..
:param resolution --- 'highest', 'lowest', 'average', 'user'.
:param outputBounds --- output bounds as (minX, minY, maxX, maxY) in target SRS.
:param xRes, yRes --- output resolution in target SRS.
:param targetAlignedPixels --- whether to force output bounds to be multiple of output resolution.
:param bandList --- array of band numbers (index start at 1).
:param addAlpha --- whether to add an alpha mask band to the VRT when the source raster have none.
:param resampleAlg --- resampling mode.
:param outputSRS --- assigned output SRS.
:param allowProjectionDifference --- whether to accept input datasets have not the same projection. Note: they will *not* be reprojected.
:param srcNodata --- source nodata value(s).
:param callback --- callback method.
:param callback_data --- user data for callback.
:return: gdal.DataSet(pathVRT)
"""

Benjamin Jakimow
committed
if len(self.mBands) == 0:
print('No VRT Inputs defined.')
return None

benjamin.jakimow@geo.hu-berlin.de
committed
assert os.path.splitext(pathVRT)[-1].lower() == '.vrt'
_kwds = dict()
supported = ['options','resolution','outputBounds','xRes','yRes','targetAlignedPixels','addAlpha','resampleAlg',
'outputSRS','allowProjectionDifference','srcNodata','VRTNodata','hideNodata','callback', 'callback_data']
for k in kwds.keys():
if k in supported:
_kwds[k] = kwds[k]

benjamin.jakimow@geo.hu-berlin.de
committed
if 'resampleAlg' not in _kwds:
_kwds['resampleAlg'] = resampleAlg

benjamin.jakimow@geo.hu-berlin.de
committed
if isinstance(self.mOutputBounds, RasterBounds):
bounds = self.mOutputBounds.polygon
xmin, ymin,xmax, ymax = bounds
_kwds['outputBounds'] = (xmin, ymin,xmax, ymax)

benjamin.jakimow@geo.hu-berlin.de
committed
dirVrt = os.path.dirname(pathVRT)
dirWarpedVRT = os.path.join(dirVrt, 'WarpedVRTs')
if not os.path.isdir(dirVrt):
os.mkdir(dirVrt)
srcLookup = dict()

benjamin.jakimow@geo.hu-berlin.de
committed
for i, pathSrc in enumerate(self.sourceRaster()):
dsSrc = gdal.Open(pathSrc)
assert isinstance(dsSrc, gdal.Dataset)
band = dsSrc.GetRasterBand(1)
noData = band.GetNoDataValue()
if noData and srcNodata is None:
srcNodata = noData

benjamin.jakimow@geo.hu-berlin.de
committed
crs = QgsCoordinateReferenceSystem(dsSrc.GetProjection())
if crs == self.mCrs:
srcLookup[pathSrc] = pathSrc
else:
if not os.path.isdir(dirWarpedVRT):
os.mkdir(dirWarpedVRT)
pathVRT2 = os.path.join(dirWarpedVRT, 'warped.{}.vrt'.format(os.path.basename(pathSrc)))
wops = gdal.WarpOptions(format='VRT',
dstSRS=self.mCrs.toWkt())
tmp = gdal.Warp(pathVRT2, dsSrc, options=wops)
assert isinstance(tmp, gdal.Dataset)
tmp = None
srcLookup[pathSrc] = pathVRT2
srcFiles = [srcLookup[src] for src in self.sourceRaster()]
#1. build a temporary VRT that described the spatial shifts of all input sources
gdal.BuildVRT(pathVRT, srcFiles, options=vro)
dsVRTDst = gdal.Open(pathVRT)
assert isinstance(dsVRTDst, gdal.Dataset)

benjamin.jakimow@geo.hu-berlin.de
committed
assert len(srcLookup) == dsVRTDst.RasterCount
ns, nl = dsVRTDst.RasterXSize, dsVRTDst.RasterYSize
gt = dsVRTDst.GetGeoTransform()
crs = dsVRTDst.GetProjectionRef()
eType = dsVRTDst.GetRasterBand(1).DataType
SOURCE_TEMPLATES = dict()
for i, srcFile in enumerate(srcFiles):
vrt_sources = dsVRTDst.GetRasterBand(i+1).GetMetadata('vrt_sources')
assert len(vrt_sources) == 1
srcXML = vrt_sources.values()[0]
assert os.path.basename(srcFile)+'</SourceFilename>' in srcXML
assert '<SourceBand>1</SourceBand>' in srcXML
SOURCE_TEMPLATES[srcFile] = srcXML
dsVRTDst = None
#remove the temporary VRT, we don't need it any more
os.remove(pathVRT)
#2. build final VRT from scratch
drvVRT = gdal.GetDriverByName('VRT')
assert isinstance(drvVRT, gdal.Driver)
dsVRTDst = drvVRT.Create(pathVRT, ns, nl,0, eType=eType)
#2.1. set general properties
assert isinstance(dsVRTDst, gdal.Dataset)
dsVRTDst.SetProjection(crs)
dsVRTDst.SetGeoTransform(gt)
#2.2. add virtual bands

benjamin.jakimow@geo.hu-berlin.de
committed
for i, vBand in enumerate(self.mBands):
assert isinstance(vBand, VRTRasterBand)
assert dsVRTDst.AddBand(eType, options=['subClass=VRTSourcedRasterBand']) == 0
vrtBandDst = dsVRTDst.GetRasterBand(i+1)
assert isinstance(vrtBandDst, gdal.Band)

benjamin.jakimow@geo.hu-berlin.de
committed
vrtBandDst.SetDescription(str(vBand.name()))
md = {}
#add all input sources for this virtual band
for iSrc, sourceInfo in enumerate(vBand.sources):

benjamin.jakimow@geo.hu-berlin.de
committed
assert isinstance(sourceInfo, VRTRasterInputSourceBand)
bandIndex = sourceInfo.mBandIndex
xml = SOURCE_TEMPLATES[srcLookup[sourceInfo.mPath]]
xml = re.sub('<SourceBand>1</SourceBand>','<SourceBand>{}</SourceBand>'.format(bandIndex+1), xml)
md['source_{}'.format(iSrc)] = xml
vrtBandDst.SetMetadata(md,'vrt_sources')
vrtBandDst.ComputeBandStats(1)
dsVRTDst = None
#check if we get what we like to get
dsCheck = gdal.Open(pathVRT)
s = ""
def __repr__(self):
info = ['VirtualRasterBuilder: {} bands, {} source files'.format(

benjamin.jakimow@geo.hu-berlin.de
committed
len(self.mBands), len(self.sourceRaster()))]
for vBand in self.mBands:
info.append(str(vBand))
return '\n'.join(info)

benjamin.jakimow@geo.hu-berlin.de
committed

benjamin.jakimow@geo.hu-berlin.de
committed
return len(self.mBands)

benjamin.jakimow@geo.hu-berlin.de
committed
return self.mBands[slice]
def __delitem__(self, slice):
self.removeVirtualBands(self[slice])
def __contains__(self, item):

benjamin.jakimow@geo.hu-berlin.de
committed
return item in self.mBands
def __iter__(self):
return iter(self.mClasses)

benjamin.jakimow@geo.hu-berlin.de
committed
class VRTRasterVectorLayer(QgsVectorLayer):
def __init__(self, vrtRaster, crs=None):
assert isinstance(vrtRaster, VRTRaster)
if crs is None:
crs = QgsCoordinateReferenceSystem('EPSG:4326')
uri = 'polygon?crs={}'.format(crs.authid())

benjamin.jakimow@geo.hu-berlin.de
committed
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
super(VRTRasterVectorLayer, self).__init__(uri, 'VRTRaster', 'memory', False)
self.mCrs = crs
self.mVRTRaster = vrtRaster
#initialize fields
assert self.startEditing()
# standard field names, types, etc.
fieldDefs = [('oid', QVariant.Int, 'integer'),
('type', QVariant.String, 'string'),
('name', QVariant.String, 'string'),
('path', QVariant.String, 'string'),
]
# initialize fields
for fieldDef in fieldDefs:
field = QgsField(fieldDef[0], fieldDef[1], fieldDef[2])
self.addAttribute(field)
self.commitChanges()
symbol = QgsFillSymbolV2.createSimple({'style': 'no', 'color': 'red', 'outline_color':'black'})
self.rendererV2().setSymbol(symbol)
self.label().setFields(self.fields())
self.label().setLabelField(3,3)
self.mVRTRaster.sigSourceRasterAdded.connect(self.onRasterInserted)
self.mVRTRaster.sigSourceRasterRemoved.connect(self.onRasterRemoved)
self.onRasterInserted(self.mVRTRaster.sourceRaster())

benjamin.jakimow@geo.hu-berlin.de
committed
def path2feature(self, path):
for f in self.dataProvider().getFeatures():
if str(f.attribute('path')) == str(path):
return f
return None
def path2fid(self, path):
for f in self.dataProvider().getFeatures():
if str(f.attribute('path')) == str(path):
return f.id()
return None
def fid2path(self, fid):
for f in self.dataProvider().getFeatures():
if f.fid() == fid:
return f
return None

benjamin.jakimow@geo.hu-berlin.de
committed
def onRasterInserted(self, listOfNewFiles):
assert isinstance(listOfNewFiles, list)
if len(listOfNewFiles) == 0:
return

benjamin.jakimow@geo.hu-berlin.de
committed
for f in listOfNewFiles:
bounds = self.mVRTRaster.sourceRasterBounds()[f]
assert isinstance(bounds, RasterBounds)
oid = str(id(bounds))
geometry =QgsPolygonV2(bounds.polygon)
#geometry = QgsCircularStringV2(bounds.curve)

benjamin.jakimow@geo.hu-berlin.de
committed
trans = QgsCoordinateTransform(bounds.crs, self.crs())
geometry.transform(trans)
feature = QgsFeature(self.pendingFields())
#feature.setGeometry(QgsGeometry(geometry))
feature.setGeometry(QgsGeometry.fromWkt(geometry.asWkt()))
#feature.setFeatureId(int(oid))
feature.setAttribute('oid', oid)
feature.setAttribute('type', 'source file')
feature.setAttribute('name', str(os.path.basename(f)))
feature.setAttribute('path', str(f))
#feature.setValid(True)

benjamin.jakimow@geo.hu-berlin.de
committed
assert self.dataProvider().addFeatures([feature])

benjamin.jakimow@geo.hu-berlin.de
committed
self.featureAdded.emit(feature.id())

benjamin.jakimow@geo.hu-berlin.de
committed

benjamin.jakimow@geo.hu-berlin.de
committed
self.updateExtents()
assert self.commitChanges()
self.dataChanged.emit()

benjamin.jakimow@geo.hu-berlin.de
committed
def onRasterRemoved(self, files):
self.startEditing()
self.selectAll()
toRemove = []
for f in self.selectedFeatures():
if f.attribute('path') in files:
toRemove.append(f.id())
self.setSelectedFeatures(toRemove)
self.deleteSelectedFeatures()
self.commitChanges()

benjamin.jakimow@geo.hu-berlin.de
committed
self.dataChanged.emit()

benjamin.jakimow@geo.hu-berlin.de
committed

benjamin.jakimow@geo.hu-berlin.de
committed
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
def createVirtualBandMosaic(bandFiles, pathVRT):
drv = gdal.GetDriverByName('VRT')
refPath = bandFiles[0]
refDS = gdal.Open(refPath)
ns, nl, nb = refDS.RasterXSize, refDS.RasterYSize, refDS.RasterCount
noData = refDS.GetRasterBand(1).GetNoDataValue()
vrtOptions = gdal.BuildVRTOptions(
# here we can use the options known from http://www.gdal.org/gdalbuildvrt.html
separate=False
)
if len(bandFiles) > 1:
s =""
vrtDS = gdal.BuildVRT(pathVRT, bandFiles, options=vrtOptions)
vrtDS.FlushCache()
assert vrtDS.RasterCount == nb
return vrtDS
def createVirtualBandStack(bandFiles, pathVRT):
nb = len(bandFiles)
drv = gdal.GetDriverByName('VRT')
refPath = bandFiles[0]
refDS = gdal.Open(refPath)
ns, nl = refDS.RasterXSize, refDS.RasterYSize
noData = refDS.GetRasterBand(1).GetNoDataValue()
vrtOptions = gdal.BuildVRTOptions(
# here we can use the options known from http://www.gdal.org/gdalbuildvrt.html

benjamin.jakimow@geo.hu-berlin.de
committed
)
vrtDS = gdal.BuildVRT(pathVRT, bandFiles, options=vrtOptions)
vrtDS.FlushCache()
assert vrtDS.RasterCount == nb
#copy band metadata from
for i in range(nb):
band = vrtDS.GetRasterBand(i+1)
band.SetDescription(bandFiles[i])

benjamin.jakimow@geo.hu-berlin.de
committed
if noData:
band.SetNoDataValue(noData)

benjamin.jakimow@geo.hu-berlin.de
committed
return vrtDS
from timeseriesviewer.utils import loadUi
class TreeNode(QObject):

benjamin.jakimow@geo.hu-berlin.de
committed
sigWillAddChildren = pyqtSignal(QObject, int, int)
sigAddedChildren = pyqtSignal(QObject, int, int)
sigWillRemoveChildren = pyqtSignal(QObject, int, int)
sigRemovedChildren = pyqtSignal(QObject, int, int)
sigUpdated = pyqtSignal(QObject)

benjamin.jakimow@geo.hu-berlin.de
committed
def __init__(self, parentNode, name=None):
super(TreeNode, self).__init__()
self.mParent = parentNode

benjamin.jakimow@geo.hu-berlin.de
committed

benjamin.jakimow@geo.hu-berlin.de
committed
self.mName = name
self.mValues = []
self.mIcon = None
self.mToolTip = None

benjamin.jakimow@geo.hu-berlin.de
committed
if isinstance(parentNode, TreeNode):
parentNode.appendChildNodes(self)

Benjamin Jakimow
committed
def nodeIndex(self):
return self.mParent.mChildren.index(self)
def next(self):
i = self.nodeIndex()
if i < len(self.mChildren.mChildren):
return self.mParent.mChildren[i+1]
else:
return None
def previous(self):
i = self.nodeIndex()
if i > 0:
return self.mParent.mChildren[i - 1]
else:
return None

benjamin.jakimow@geo.hu-berlin.de
committed
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
def detach(self):
"""
Detaches this TreeNode from its parent TreeNode
:return:
"""
if isinstance(self.mParent, TreeNode):
self.mParent.mChildren.remove(self)
self.setParentNode(None)
def appendChildNodes(self, listOfChildNodes):
self.insertChildNodes(len(self.mChildren), listOfChildNodes)
def insertChildNodes(self, index, listOfChildNodes):
assert index <= len(self.mChildren)
if isinstance(listOfChildNodes, TreeNode):
listOfChildNodes = [listOfChildNodes]
assert isinstance(listOfChildNodes, list)
l = len(listOfChildNodes)
idxLast = index+l-1
self.sigWillAddChildren.emit(self, index, idxLast)
for i, node in enumerate(listOfChildNodes):
assert isinstance(node, TreeNode)
node.mParent = self
# connect node signals
node.sigWillAddChildren.connect(self.sigWillAddChildren)
node.sigAddedChildren.connect(self.sigAddedChildren)
node.sigWillRemoveChildren.connect(self.sigWillRemoveChildren)
node.sigRemovedChildren.connect(self.sigRemovedChildren)
node.sigUpdated.connect(self.sigUpdated)
self.mChildren.insert(index+i, node)
self.sigAddedChildren.emit(self, index, idxLast)
def removeChildNode(self, node):
assert node in self.mChildren
i = self.mChildren.index(node)
self.removeChildNodes(i, 1)
def removeChildNodes(self, row, count):
if row < 0 or count <= 0:
return False
rowLast = row + count - 1
if rowLast >= self.childCount():
return False
self.sigWillRemoveChildren.emit(self, row, rowLast)
to_remove = self.childNodes()[row:rowLast+1]
for n in to_remove:
self.mChildren.remove(n)
#n.mParent = None
self.sigRemovedChildren.emit(self, row, rowLast)
def setToolTip(self, toolTip):
self.mToolTip = toolTip
def toolTip(self):
return self.mToolTip

benjamin.jakimow@geo.hu-berlin.de
committed
def parentNode(self):

benjamin.jakimow@geo.hu-berlin.de
committed
def setParentNode(self, treeNode):
assert isinstance(treeNode, TreeNode)
self.mParent = treeNode
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
def setIcon(self, icon):
self.mIcon = icon
def icon(self):
return self.mIcon
def setName(self, name):
self.mName = name
def name(self):
return self.mName
def contextMenu(self):
return None
def setValues(self, listOfValues):
if not isinstance(listOfValues, list):
listOfValues = [listOfValues]
self.mValues = listOfValues[:]
def values(self):
return self.mValues[:]
def childCount(self):
return len(self.mChildren)

benjamin.jakimow@geo.hu-berlin.de
committed
def childNodes(self):
return self.mChildren[:]
def findChildNodes(self, type, recursive=True):
results = []
for node in self.mChildren:
if isinstance(node, type):
results.append(node)
if recursive:
results.extend(node.findChildNodes(type, recursive=True))
return results
class SourceRasterFileNode(TreeNode):
def __init__(self, parentNode, path):
super(SourceRasterFileNode, self).__init__(parentNode)
self.mPath = path
self.setName(os.path.basename(path))

benjamin.jakimow@geo.hu-berlin.de
committed
srcNode = TreeNode(self, name='Path')
srcNode.setValues(path)
#populate metainfo
ds = gdal.Open(path)
assert isinstance(ds, gdal.Dataset)

benjamin.jakimow@geo.hu-berlin.de
committed

benjamin.jakimow@geo.hu-berlin.de
committed

benjamin.jakimow@geo.hu-berlin.de
committed
crsNode = TreeNode(self, name='CRS')
crsNode.setIcon(QIcon(':/timeseriesviewer/icons/CRS.png'))

benjamin.jakimow@geo.hu-berlin.de
committed
crs = osr.SpatialReference()
crs.ImportFromWkt(ds.GetProjection())
authInfo = '{}:{}'.format(crs.GetAttrValue('AUTHORITY',0), crs.GetAttrValue('AUTHORITY',1))
crsNode.setValues([authInfo,crs.ExportToWkt()])
self.bandNode = TreeNode(None, name='Bands')
for b in range(ds.RasterCount):
band = ds.GetRasterBand(b+1)

benjamin.jakimow@geo.hu-berlin.de
committed
inputSource = VRTRasterInputSourceBand(path, b)
inputSource.mBandName = band.GetDescription()
if inputSource.mBandName in [None,'']:
inputSource.mBandName = '{}'.format(b + 1)
inputSource.mNoData = band.GetNoDataValue()
SourceRasterBandNode(self.bandNode, inputSource)
self.bandNode.setParentNode(self)
self.appendChildNodes(self.bandNode)
def sourceBands(self):
return [n.mSrcBand for n in self.bandNode.mChildren if isinstance(n, SourceRasterBandNode)]
class SourceRasterBandNode(TreeNode):

benjamin.jakimow@geo.hu-berlin.de
committed
def __init__(self, parentNode, vrtRasterInputSourceBand):
assert isinstance(vrtRasterInputSourceBand, VRTRasterInputSourceBand)
super(SourceRasterBandNode, self).__init__(parentNode)
self.setIcon(QIcon(":/timeseriesviewer/icons/mIconRaster.png"))

benjamin.jakimow@geo.hu-berlin.de
committed
self.mSrcBand = vrtRasterInputSourceBand
self.setName(self.mSrcBand.mBandName)
#self.setValues([self.mSrcBand.mPath])
self.setToolTip('band {}:{}'.format(self.mSrcBand.mBandIndex+1, self.mSrcBand.mPath))
class VRTRasterNode(TreeNode):
def __init__(self, parentNode, vrtRaster):
assert isinstance(vrtRaster, VRTRaster)
super(VRTRasterNode, self).__init__(parentNode)
self.mVRTRaster = vrtRaster
self.mVRTRaster.sigBandInserted.connect(self.onBandInserted)
self.mVRTRaster.sigBandRemoved.connect(self.onBandRemoved)

benjamin.jakimow@geo.hu-berlin.de
committed
def onBandInserted(self, index, vrtRasterBand):
assert isinstance(vrtRasterBand, VRTRasterBand)
i = vrtRasterBand.bandIndex()
assert i == index
node = VRTRasterBandNode(None, vrtRasterBand)
self.insertChildNodes(i, [node])
def onBandRemoved(self, removedIdx):
self.removeChildNodes(removedIdx, 1)
class VRTRasterBandNode(TreeNode):
def __init__(self, parentNode, virtualBand):

benjamin.jakimow@geo.hu-berlin.de
committed
assert isinstance(virtualBand, VRTRasterBand)
super(VRTRasterBandNode, self).__init__(parentNode)
self.mVirtualBand = virtualBand
self.setName(virtualBand.name())
self.setIcon(QIcon(":/timeseriesviewer/icons/mIconVirtualRaster.png"))

benjamin.jakimow@geo.hu-berlin.de
committed
#self.nodeBands = TreeNode(self, name='Input Bands')
#self.nodeBands.setToolTip('Source bands contributing to this virtual raster band')
self.nodeBands = self
virtualBand.sigNameChanged.connect(self.setName)

benjamin.jakimow@geo.hu-berlin.de
committed
virtualBand.sigSourceInserted.connect(lambda _, src: self.onSourceInserted(src))
virtualBand.sigSourceRemoved.connect(self.onSourceRemoved)
for src in self.mVirtualBand.sources:
self.onSourceInserted(src)

benjamin.jakimow@geo.hu-berlin.de
committed
def onSourceInserted(self, inputSource):
assert isinstance(inputSource, VRTRasterInputSourceBand)
assert inputSource.virtualBand() == self.mVirtualBand
i = self.mVirtualBand.sources.index(inputSource)

benjamin.jakimow@geo.hu-berlin.de
committed
node = VRTRasterInputSourceBandNode(None, inputSource)
self.nodeBands.insertChildNodes(i, node)

benjamin.jakimow@geo.hu-berlin.de
committed
def onSourceRemoved(self, row, inputSource):
assert isinstance(inputSource, VRTRasterInputSourceBand)
node = self.nodeBands.childNodes()[row]
if node.mSrc != inputSource:
s = ""
self.nodeBands.removeChildNode(node)

benjamin.jakimow@geo.hu-berlin.de
committed
class VRTRasterInputSourceBandNode(TreeNode):
def __init__(self, parentNode, vrtRasterInputSourceBand):
assert isinstance(vrtRasterInputSourceBand, VRTRasterInputSourceBand)
super(VRTRasterInputSourceBandNode, self).__init__(parentNode)
self.setIcon(QIcon(":/timeseriesviewer/icons/mIconRaster.png"))

benjamin.jakimow@geo.hu-berlin.de
committed
self.mSrc = vrtRasterInputSourceBand
name = '{}:{}'.format(self.mSrc.mBandIndex+1, os.path.basename(self.mSrc.mPath))
self.setName(name)
#self.setValues([self.mSrc.mPath, self.mSrc.mBandIndex])
def sourceBand(self):
return self.mSrc
class TreeView(QTreeView):
def __init__(self, *args, **kwds):
super(TreeView, self).__init__(*args, **kwds)
class TreeModel(QAbstractItemModel):
def __init__(self, parent=None, rootNode = None):
super(TreeModel, self).__init__(parent)
self.mColumnNames = ['Node','Value']

benjamin.jakimow@geo.hu-berlin.de
committed
self.mRootNode = rootNode if isinstance(rootNode, TreeNode) else TreeNode(None)
self.mRootNode.sigWillAddChildren.connect(self.nodeWillAddChildren)
self.mRootNode.sigAddedChildren.connect(self.nodeAddedChildren)
self.mRootNode.sigWillRemoveChildren.connect(self.nodeWillRemoveChildren)
self.mRootNode.sigRemovedChildren.connect(self.nodeRemovedChildren)
self.mRootNode.sigUpdated.connect(self.nodeUpdated)
self.mTreeView = None
if isinstance(parent, QTreeView):
self.connectTreeView(parent)
def nodeWillAddChildren(self, node, idx1, idxL):
idxNode = self.node2idx(node)
self.beginInsertRows(idxNode, idx1, idxL)
def nodeAddedChildren(self, node, idx1, idxL):
self.endInsertRows()
#for i in range(idx1, idxL+1):
for n in node.childNodes():
self.setColumnSpan(node)