From 30d602bdaf8d10b2046bd52e36da9e064d967f89 Mon Sep 17 00:00:00 2001
From: Benjamin Jakimow <no.email>
Date: Fri, 9 Dec 2016 21:27:54 +0100
Subject: [PATCH] refactoring
uses QgsRasterLayers
---
timeseriesviewer/timeseries.py | 763 +++++++++++++++++++++++++++++++++
1 file changed, 763 insertions(+)
create mode 100644 timeseriesviewer/timeseries.py
diff --git a/timeseriesviewer/timeseries.py b/timeseriesviewer/timeseries.py
new file mode 100644
index 00000000..ff640acd
--- /dev/null
+++ b/timeseriesviewer/timeseries.py
@@ -0,0 +1,763 @@
+from __future__ import absolute_import
+import six, sys, os, gc, re, collections, site, inspect, time, traceback, copy
+import bisect
+from osgeo import gdal, ogr
+
+from qgis import *
+from qgis.core import *
+from qgis.gui import *
+from PyQt4.QtGui import *
+from PyQt4.QtCore import *
+
+import numpy as np
+
+
+def transformGeometry(geom, crsSrc, crsDst, trans=None):
+ if trans is None:
+ assert isinstance(crsSrc, QgsCoordinateReferenceSystem)
+ assert isinstance(crsDst, QgsCoordinateReferenceSystem)
+ return transformGeometry(geom, None, None, trans=QgsCoordinateTransform(crsSrc, crsDst))
+ else:
+ assert isinstance(trans, QgsCoordinateTransform)
+ return trans.transform(geom)
+
+
+
+
+class SensorInstrument(object):
+
+ INSTRUMENTS = dict()
+ INSTRUMENTS = {(6, 30., 30.): 'Landsat Legacy' \
+ , (7, 30., 30.): 'L8 OLI' \
+ , (4, 10., 10.): 'S2 MSI 10m' \
+ , (6, 20., 20.): 'S2 MSI 20m' \
+ , (3, 30., 30.): 'S2 MSI 60m' \
+ , (3, 30., 30.): 'S2 MSI 60m' \
+ , (5, 5., 5.): 'RE 5m' \
+ }
+
+ @staticmethod
+ def fromRasterLayer(lyr):
+
+ assert isinstance(lyr, QgsRasterLayer)
+ nb = lyr.bandCount()
+ sx = lyr.rasterUnitsPerPixelX()
+ sy = lyr.rasterUnitsPerPixelY()
+
+ bandNames = [lyr.bandName(i) for i in range(1, nb+1)]
+
+ return SensorInstrument(nb, sx, sy,
+ band_names=bandNames,
+ wavelengths=None,
+ sensor_name=None)
+
+ def fromGDALDataSet(self, ds):
+ assert isinstance(ds, gdal.Dataset)
+ nb = ds.RasterCount
+
+
+ """
+ Describes a Sensor Configuration
+ """
+ def __init__(self,nb, px_size_x,px_size_y, band_names=None, wavelengths=None, sensor_name=None):
+ assert nb >= 1
+
+ self.TS = None
+ self.nb = nb
+ self.px_size_x = float(abs(px_size_x))
+ self.px_size_y = float(abs(px_size_y))
+
+ assert self.px_size_x > 0
+ assert self.px_size_y > 0
+
+ if band_names is not None:
+ assert len(band_names) == nb
+ else:
+ band_names = ['Band {}'.format(b+1) for b in range(nb)]
+
+ self.bandNames = band_names
+
+ if wavelengths is not None:
+ assert len(wavelengths) == nb
+
+ self.wavelengths = wavelengths
+
+ if sensor_name is None:
+ id = (self.nb, self.px_size_x, self.px_size_y)
+ sensor_name = SensorInstrument.INSTRUMENTS.get(
+ id,
+ '{}band@{}m'.format(self.nb, self.px_size_x))
+
+ self.sensorName = sensor_name
+
+ self.hashvalue = hash(','.join(self.bandNames))
+
+ def __eq__(self, other):
+ return self.nb == other.nb and \
+ self.px_size_x == other.px_size_x and \
+ self.px_size_y == other.px_size_y
+
+ def __hash__(self):
+ return self.hashvalue
+
+ def __repr__(self):
+ return self.sensorName
+
+ def getDescription(self):
+ info = []
+ info.append(self.sensorName)
+ info.append('{} Bands'.format(self.nb))
+ info.append('Band\tName\tWavelength')
+ for b in range(self.nb):
+ if self.wavelengths:
+ wl = str(self.wavelengths[b])
+ else:
+ wl = 'unknown'
+ info.append('{}\t{}\t{}'.format(b + 1, self.bandNames[b], wl))
+
+ return '\n'.join(info)
+
+
+class TSDLoaderSignals(QObject):
+
+ sigRasterLayerLoaded = pyqtSignal(QgsRasterLayer)
+ sigFinished = pyqtSignal()
+
+class TSDLoader(QRunnable):
+ """
+ Runnable to load QgsRasterLayers from a parallel thread
+ """
+ def __init__(self, tsd_paths):
+ super(TSDLoader, self).__init__()
+ self.signals = TSDLoaderSignals()
+
+ self.paths = list()
+ for p in tsd_paths:
+ if not (isinstance(p, tuple) or isinstance(p, list)):
+ p = [p]
+ self.paths.append(p)
+
+
+ def run(self):
+ lyrs = []
+ for path in self.paths:
+ TSD = TimeSeriesDatum(path)
+ lyr = QgsRasterLayer(path)
+ if lyr:
+ lyrs.append(lyr)
+ self.signals.sigRasterLayerLoaded.emit(lyr)
+ dprint('{} loaded'.format(path))
+ else:
+ dprint('Failed to load {}'.format(path))
+ self.signals.sigFinished.emit()
+ #return lyrs
+
+
+class TimeSeries(QObject):
+ datumAdded = pyqtSignal(name='datumAdded')
+
+ #fire when a new sensor configuration is added
+ sensorAdded = pyqtSignal(SensorInstrument, name='sensorAdded')
+
+ changed = pyqtSignal()
+
+ progress = pyqtSignal(int,int,int, name='progress')
+ chipLoaded = pyqtSignal(object, name='chiploaded')
+ closed = pyqtSignal()
+ error = pyqtSignal(object)
+
+ def __init__(self, imageFiles=None, maskFiles=None):
+ QObject.__init__(self)
+
+ #define signals
+
+ #fire when a new TSD is added
+
+
+ #self.data = collections.OrderedDict()
+ self.data = list()
+
+ self.CHIP_BUFFER=dict()
+
+ self.shape = None
+
+ self.Sensors = collections.OrderedDict()
+
+ self.Pool = None
+
+ if imageFiles is not None:
+ self.addFiles(imageFiles)
+ if maskFiles is not None:
+ self.addMasks(maskFiles)
+
+ _sep = ';'
+
+
+ def loadFromFile(self, path):
+
+ images = []
+ masks = []
+ with open(path, 'r') as f:
+ lines = f.readlines()
+ for l in lines:
+ if re.match('^[ ]*[;#&]', l):
+ continue
+
+ parts = re.split('[\n'+TimeSeries._sep+']', l)
+ parts = [p for p in parts if p != '']
+ images.append(parts[0])
+ if len(parts) > 1:
+ masks.append(parts[1])
+
+ self.addFiles(images)
+ self.addMasks(masks)
+
+
+ def saveToFile(self, path):
+ if path is None or len(path) == 0:
+ return
+
+ lines = []
+ lines.append('#Time series definition file: {}'.format(np.datetime64('now').astype(str)))
+ lines.append('#<image path>[;<mask path>]')
+ for TSD in self.data.values():
+
+ line = TSD.pathImg
+ if TSD.pathMsk is not None:
+ line += TimeSeries._sep + TSD.pathMsk
+
+ lines.append(line)
+
+ lines = [l+'\n' for l in lines]
+
+ print('Write {}'.format(path))
+ with open(path, 'w') as f:
+ f.writelines(lines)
+
+
+ def getSRS(self):
+ if len(self.data) == 0:
+ return 0
+ else:
+ TSD = self.data[0]
+ return TSD.getSpatialReference()
+
+ def getWKT(self):
+ srs = self.getSRS()
+ return srs.ExportToWkt()
+
+ def getSceneCenter(self, srs=None):
+
+ if srs is None:
+ srs = self.getSRS()
+
+ bbs = list()
+ for S, TSDs in self.Sensors.items():
+ x = []
+ y = []
+ for TSD in TSDs:
+ bb = TSD.getBoundingBox(srs)
+ x.extend([c[0] for c in bb])
+ y.extend([c[1] for c in bb])
+
+ return None
+ pass
+
+ def getMaxExtent(self, srs=None):
+
+
+ if srs is None:
+ srs = self.getSRS()
+ bb = None
+ for TSD in self.data:
+ if bb is None:
+ bb = TSD.getBoundingBox(srs)
+ else:
+ bb.unionRect(TSD.getBoundingBox(srs))
+
+ return bb
+
+ def getObservationDates(self):
+ return [tsd.getDate() for tsd in self.data]
+
+ def getTSDs(self, date_of_interest=None):
+ if date_of_interest:
+ tsds = [tsd for tsd in self.data if tsd.getDate() == date_of_interest]
+ else:
+ tsds = self.data
+ return tsds
+
+ def _callback_error(self, error):
+ six.print_(error, file=sys.stderr)
+ self.error.emit(error)
+ self._callback_progress()
+
+ def _callback_spatialchips(self, results):
+ self.chipLoaded.emit(results)
+ self._callback_progress()
+
+ def _callback_progress(self):
+ self._callback_progress_done += 1
+ self.progress.emit(0, self._callback_progress_done, self._callback_progress_max)
+
+ if self._callback_progress_done >= self._callback_progress_max:
+ self._callback_progress_done = 0
+ self._callback_progress_max = 0
+ self.progress.emit(0,0,1)
+
+ def getSpatialChips_parallel(self, bbWkt, srsWkt, TSD_band_list, ncpu=1):
+ assert type(bbWkt) is str and type(srsWkt) is str
+
+ import multiprocessing
+ if self.Pool is not None:
+ self.Pool.terminate()
+
+ self.Pool = multiprocessing.Pool(processes=ncpu)
+
+
+ self._callback_progress_max = len(TSD_band_list)
+ self._callback_progress_done = 0
+
+ for T in TSD_band_list:
+
+ TSD = copy.deepcopy(T[0])
+ bands = T[1]
+ #TSD = pickle.dumps(self.data[date])
+ args = (TSD, bbWkt, srsWkt)
+ kwds = {'bands':bands}
+
+ if six.PY3:
+ self.Pool.apply_async(PFunc_TimeSeries_getSpatialChip, \
+ args=args, kwds=kwds, \
+ callback=self._callback_spatialchips, error_callback=self._callback_error)
+ else:
+ self.Pool.apply_async(PFunc_TimeSeries_getSpatialChip, \
+ args, kwds, self._callback_spatialchips)
+
+ s = ""
+
+ pass
+
+ def getImageChips(self, xy, size=50, bands=[4,5,6], dates=None):
+ chipCollection = collections.OrderedDict()
+
+ if dates is None:
+ dates = self.data.keys()
+ for date in dates:
+ TSD = self.data[date]
+ chipCollection[date] = TSD.readImageChip(xy, size=size, bands=bands)
+
+ return chipCollection
+
+ def addMasks(self, files, raise_errors=True, mask_value=0, exclude_mask_value=True):
+ assert isinstance(files, list)
+ l = len(files)
+
+ self.progress.emit(0,0,l)
+ for i, file in enumerate(files):
+
+ try:
+ self.addMask(file, raise_errors=raise_errors, mask_value=mask_value, exclude_mask_value=exclude_mask_value, _quiet=True)
+ except:
+ pass
+
+ self.progress.emit(0,i+1,l)
+
+ self.progress.emit(0,0,1)
+ self.changed.emit()
+
+ def addMask(self, pathMsk, raise_errors=True, mask_value=0, exclude_mask_value=True, _quiet=False):
+ print('Add mask {}...'.format(pathMsk))
+ ds = getDS(pathMsk)
+ date = getImageDate(ds)
+
+ if date in self.data.keys():
+ TSD = self.data[date]
+
+ if not _quiet:
+ self.changed.emit()
+
+ return TSD.setMask(pathMsk, raise_errors=raise_errors, mask_value=mask_value, exclude_mask_value=exclude_mask_value)
+ else:
+ info = 'TimeSeries does not contain date {} {}'.format(date, pathMsk)
+ if raise_errors:
+ raise Exception(info)
+ else:
+ six.print_(info, file=sys.stderr)
+ return False
+
+ def removeAll(self):
+ self.clear()
+
+ def clear(self):
+ self.Sensors.clear()
+ del self.data[:]
+ self.changed.emit()
+
+
+ def removeDates(self, TSDs):
+ for TSD in TSDs:
+ self.removeTSD(TSD, _quiet=True)
+ self.changed.emit()
+
+ def removeTSD(self, TSD, _quiet=False):
+
+ assert type(TSD) is TimeSeriesDatum
+ S = TSD.sensor
+ self.Sensors[S].remove(TSD)
+ self.data.pop(TSD, None)
+ if len(self.Sensors[S]) == 0:
+ self.Sensors.pop(S)
+ if not _quiet:
+ self.changed.emit()
+
+
+
+ def addFile(self, pathImg, pathMsk=None, _quiet=False):
+ six.print_(pathImg)
+ six.print_('Add image {}...'.format(pathImg))
+ TSD = TimeSeriesDatum(pathImg, pathMsk=pathMsk)
+ self.addTimeSeriesDatum(TSD)
+
+ def addTimeSeriesDatum(self, TSD):
+
+ try:
+ sensorAdded = False
+ existingSensors = list(self.Sensors.keys())
+ if TSD.sensor not in existingSensors:
+ self.Sensors[TSD.sensor] = list()
+ sensorAdded = True
+ else:
+ TSD.sensor = existingSensors[existingSensors.index(TSD.sensor)]
+
+ if TSD in self.data:
+ six.print_('Time series datum already added: {}'.format(str(TSD)), file=sys.stderr)
+ else:
+ self.Sensors[TSD.sensor].append(TSD)
+
+ #insert sorted
+ bisect.insort(self.data, TSD)
+ #self.data[TSD] = TSD
+
+ if sensorAdded:
+ self.sensorAdded.emit(TSD.sensor)
+
+
+ except:
+ exc_type, exc_value, exc_traceback = sys.exc_info()
+ traceback.print_exception(exc_type, exc_value, exc_traceback, limit=2)
+ six.print_('Unable to add {}'.format(file), file=sys.stderr)
+ pass
+
+
+
+ def addFilesAsync(self, files):
+ assert isinstance(files, list)
+
+
+ def addFiles(self, files):
+ assert isinstance(files, list)
+ l = len(files)
+ assert l > 0
+
+ self.progress.emit(0,0,l)
+ for i, file in enumerate(files):
+ self.addFile(file, _quiet=True)
+ self.progress.emit(0,i+1,l)
+
+ self.progress.emit(0,0,1)
+ self.datumAdded.emit()
+ self.changed.emit()
+
+
+
+
+
+ def __len__(self):
+ return len(self.data)
+
+ def __repr__(self):
+ info = []
+ info.append('TimeSeries:')
+ l = len(self)
+ info.append(' Scenes: {}'.format(l))
+
+ if l > 0:
+ keys = list(self.data.keys())
+ info.append(' Range: {} to {}'.format(keys[0], keys[-1]))
+ return '\n'.join(info)
+
+regAcqDate = re.compile(r'acquisition[ ]*(time|date|day)', re.I)
+regLandsatSceneID = re.compile(r"L[EMCT][1234578]{1}[12]\d{12}[a-zA-Z]{3}\d{2}")
+
+def getImageDate2(lyr):
+ assert isinstance(lyr, QgsRasterLayer)
+ mdLines = str(lyr.metadata()).splitlines()
+ date = None
+ #find date in metadata
+ for line in [l for l in mdLines if regAcqDate.search(l)]:
+ date = parseAcquisitionDate(line)
+ if date:
+ return date
+ #find date in filename
+ dn, fn = os.path.split(str(lyr.dataProvider().dataSourceUri()))
+ date = parseAcquisitionDate(fn)
+ if date: return date
+
+ #find date in file directory path
+ date = parseAcquisitionDate(date)
+
+ return date
+
+class TimeSeriesDatum(QObject):
+ """
+ Collects all data sets related to one sensor
+ """
+
+ def __init__(self, pathImg, pathMsk=None):
+ self.pathImg = pathImg
+ self.pathMsk = None
+
+ self.lyrImg = QgsRasterLayer(pathImg, os.path.basename(pathImg), False)
+ self.crs = self.lyrImg.dataProvider().crs()
+ self.sensor = SensorInstrument.fromRasterLayer(self.lyrImg)
+
+ self.date = getImageDate2(self.lyrImg)
+ assert self.date is not None, 'Unable to find acquisition date of {}'.format(pathImg)
+
+ self.ns = self.lyrImg.width()
+ self.nl = self.lyrImg.height()
+ self.nb = self.lyrImg.bandCount()
+ self.srs_wkt = str(self.crs.toWkt())
+
+
+ if pathMsk:
+ self.setMask(pathMsk)
+
+ def getdtype(self):
+ return gdal_array.GDALTypeCodeToNumericTypeCode(self.etype)
+
+ def getDate(self):
+ return np.datetime64(self.date)
+
+
+ def getSpatialReference(self):
+ return self.crs
+ srs = osr.SpatialReference()
+ srs.ImportFromWkt(self.srs_wkt)
+ return srs
+
+
+ def getBoundingBox(self, srs=None):
+
+ bbox = self.lyrImg.extent()
+ if srs:
+ assert isinstance(srs, QgsCoordinateReferenceSystem)
+ bbox = transformGeometry(bbox, self.crs, srs)
+ return bbox
+
+
+ def setMask(self, pathMsk, raise_errors=True, mask_value=0, exclude_mask_value=True):
+ dsMsk = gdal.Open(pathMsk)
+ mskDate = getImageDate(dsMsk)
+
+
+ errors = list()
+ if mskDate and mskDate != self.getDate():
+ errors.append('Mask date differs from image date')
+ if self.ns != dsMsk.RasterXSize or self.nl != dsMsk.RasterYSize:
+ errors.append('Spatial size differs')
+ if dsMsk.RasterCount != 1:
+ errors.append('Mask has > 1 bands')
+
+ projImg = self.getSpatialReference()
+ projMsk = osr.SpatialReference()
+ projMsk.ImportFromWkt(dsMsk.GetProjection())
+
+ if not projImg.IsSame(projMsk):
+ errors.append('Spatial Reference differs from image')
+ if self.gt != list(dsMsk.GetGeoTransform()):
+ errors.append('Geotransformation differs from image')
+
+ if len(errors) > 0:
+ errors.insert(0, 'pathImg:{} \npathMsk:{}'.format(self.pathImg, pathMsk))
+ errors = '\n'.join(errors)
+ if raise_errors:
+ raise Exception(errors)
+ else:
+ six.print_(errors, file=sys.stderr)
+ return False
+ else:
+ self.pathMsk = pathMsk
+ self.mask_value = mask_value
+ self.exclude_mask_value = exclude_mask_value
+
+ return True
+
+ def readSpatialChip(self, geometry, srs=None, bands=[4,5,3]):
+
+ srs_img = osr.SpatialReference()
+ srs_img.ImportFromWkt(self.srs_wkt)
+
+
+ if type(geometry) is ogr.Geometry:
+ g_bb = geometry
+ srs_bb = g_bb.GetSpatialReference()
+ else:
+ assert srs is not None and type(srs) in [str, osr.SpatialReference]
+ if type(srs) is str:
+ srs_bb = osr.SpatialReference()
+ srs_bb.ImportFromWkt(srs)
+ else:
+ srs_bb = srs.Clone()
+ g_bb = ogr.CreateGeometryFromWkt(geometry, srs_bb)
+
+ assert srs_bb is not None and g_bb is not None
+ assert g_bb.GetGeometryName() == 'POLYGON'
+
+ if not srs_img.IsSame(srs_bb):
+ g_bb = g_bb.Clone()
+ g_bb.TransformTo(srs_img)
+
+ cx0,cx1,cy0,cy1 = g_bb.GetEnvelope()
+
+ ul_px = coordinate2px(self.gt, min([cx0, cx1]), max([cy0, cy1]))
+ lr_px = coordinate2px(self.gt, max([cx0, cx1]), min([cy0, cy1]))
+ lr_px = [c+1 for c in lr_px] #+1
+
+ return self.readImageChip([ul_px[0], lr_px[0]], [ul_px[1], lr_px[1]], bands=bands)
+
+ def readImageChip(self, px_x, px_y, bands=[4,5,3]):
+
+ ds = gdal.Open(self.pathImg, gdal.GA_ReadOnly)
+
+ assert len(px_x) == 2 and px_x[0] <= px_x[1]
+ assert len(px_y) == 2 and px_y[0] <= px_y[1]
+
+ ns = px_x[1]-px_x[0]+1
+ nl = px_y[1]-px_y[0]+1
+ assert ns >= 0
+ assert nl >= 0
+
+ src_ns = ds.RasterXSize
+ src_nl = ds.RasterYSize
+
+
+ chipdata = dict()
+
+
+
+ #pixel indices in source image
+ x0 = max([0, px_x[0]])
+ y0 = max([0, px_y[0]])
+ x1 = min([src_ns, px_x[1]])
+ y1 = min([src_nl, px_y[1]])
+ win_xsize = x1-x0+1
+ win_ysize = y1-y0+1
+
+ #pixel indices in image chip (ideally 0 and ns-1 or nl-1)
+ i0 = x0 - px_x[0]
+ i1 = i0 + win_xsize
+
+ j0 = y0 - px_y[0]
+ j1 = j0+ win_ysize
+
+
+
+
+ templateImg = np.zeros((nl,ns))
+ if self.nodata:
+ templateImg *= self.nodata
+
+ templateImg = templateImg.astype(self.getdtype())
+ templateMsk = np.ones((nl,ns), dtype='bool')
+
+ if win_xsize < 1 or win_ysize < 1:
+ six.print_('Selected image chip is out of raster image {}'.format(self.pathImg), file=sys.stderr)
+ for i, b in enumerate(bands):
+ chipdata[b] = np.copy(templateImg)
+
+ else:
+ for i, b in enumerate(bands):
+ band = ds.GetRasterBand(b)
+ data = np.copy(templateImg)
+ data[j0:j1,i0:i1] = band.ReadAsArray(xoff=x0, yoff=y0, win_xsize=win_xsize,win_ysize=win_ysize)
+ chipdata[b] = data
+ nodatavalue = band.GetNoDataValue()
+ if nodatavalue is not None:
+ templateMsk[j0:j1,i0:i1] = np.logical_and(templateMsk[j0:j1,i0:i1], data[j0:j1,i0:i1] != nodatavalue)
+
+ if self.pathMsk:
+ ds = gdal.Open(self.pathMsk)
+ tmp = ds.GetRasterBand(1).ReadAsArray(xoff=x0, yoff=y0, \
+ win_xsize=win_xsize,win_ysize=win_ysize) == 0
+
+ templateMsk[j0:j1,i0:i1] = np.logical_and(templateMsk[j0:j1,i0:i1], tmp)
+
+ chipdata['mask'] = templateMsk
+ return chipdata
+
+ def __repr__(self):
+ return 'TS Datum {} {}'.format(self.date, str(self.sensor))
+
+ def __cmp__(self, other):
+ return cmp(str((self.date, self.sensor)), str((other.date, other.sensor)))
+
+ def __eq__(self, other):
+ return self.date == other.date and self.sensor == other.sensor
+
+ def __hash__(self):
+ return hash((self.date,self.sensor.sensorName))
+
+
+def PFunc_TimeSeries_getSpatialChip(TSD, bbWkt, srsWkt , bands=[4,5,3]):
+
+ chipdata = TSD.readSpatialChip(bbWkt, srs=srsWkt, bands=bands)
+
+ return TSD, chipdata
+
+def px2Coordinate(gt, pxX, pxY, upper_left=True):
+ cx = gt[0] + pxX*gt[1] + pxY*gt[2]
+ cy = gt[3] + pxX*gt[4] + pxY*gt[5]
+ if not upper_left:
+ cx += gt[1]*0.5
+ cy += gt[5]*0.5
+ return cx, cy
+
+def coordinate2px(gt, cx, cy):
+ px = int((cx - gt[0]) / gt[1])
+ py = int((cy - gt[3]) / gt[5])
+ return px, py
+
+
+regYYYYDOY = re.compile(r'(19|20)\d{5}')
+regYYYYMMDD = re.compile(r'(19|20)\d{2}-\d{2}-\d{2}')
+regYYYY = re.compile(r'(19|20)\d{2}')
+
+def parseAcquisitionDate(text):
+ match = regLandsatSceneID.search(text)
+ if match:
+ id = match.group()
+ return getDateTime64FromYYYYDOY(id[9:16])
+ match = regYYYYMMDD.search(text)
+ if match:
+ return np.datetime64(match.group())
+ match = regYYYYDOY.search(text)
+ if match:
+ return getDateTime64FromYYYYDOY(match.group())
+ match = regYYYY.search(text)
+ if match:
+ return np.datetime64(match.group())
+ return None
+
+
+
+def getDateTime64FromYYYYDOY(yyyydoy):
+ return getDateTime64FromDOY(yyyydoy[0:4], yyyydoy[4:7])
+
+def getDateTime64FromDOY(year, doy):
+ if type(year) is str:
+ year = int(year)
+ if type(doy) is str:
+ doy = int(doy)
+ return np.datetime64('{:04d}-01-01'.format(year)) + np.timedelta64(doy-1, 'D')
--
GitLab