|
| 1 | +__author__ = 'Jonas Eberle <jonas.eberle@eberle-mail.de>' |
| 2 | +__version__ = '0.1' |
| 3 | + |
| 4 | +import logging |
| 5 | +from pyEOM.sources import * |
| 6 | +from pyEOM import processing |
| 7 | +from datetime import datetime |
| 8 | +import os, re, sys |
| 9 | +import subprocess |
| 10 | +import sqlite3 |
| 11 | +from pyspatialite import dbapi2 as db |
| 12 | +from datetime import datetime |
| 13 | +import shutil |
| 14 | + |
| 15 | +from pyEOM import log |
| 16 | +LOGGER = log.LOGGER |
| 17 | + |
| 18 | + |
| 19 | +class Tiles(object): |
| 20 | + database = os.path.dirname(os.path.abspath(processing.__file__))+'/data/modis_tiles.sqlite' |
| 21 | + |
| 22 | + def __init__(self): |
| 23 | + self.dbcon = db.connect(self.database) |
| 24 | + # self.dbcon.row_factory = sqlite3.Row |
| 25 | + self.dbcur = self.dbcon.cursor() |
| 26 | + |
| 27 | + def getTilesFromGeometry(self, geom): |
| 28 | + tiles = [] |
| 29 | + |
| 30 | + sel = self.dbcur.execute("SELECT 'h'||substr('00'||cast(h as int), -2,2 )||'v'||substr('00'||cast(v as int),-2,2) as tile FROM modis_tiles WHERE Intersects(GEOMETRY, GeomFromText('%s'))==1" % (geom)).fetchall() |
| 31 | + if len(sel) > 0: |
| 32 | + for tile in sel: |
| 33 | + tiles.append(tile[0]) |
| 34 | + return tiles |
| 35 | + |
| 36 | + |
| 37 | +class GEE(GEESource): |
| 38 | + |
| 39 | + def __init__(self, task, dataset): |
| 40 | + self.connect() |
| 41 | + pass |
| 42 | + |
| 43 | + def process(self): |
| 44 | + pass |
| 45 | + |
| 46 | + def ingest(self, downloadPath=None): |
| 47 | + |
| 48 | + pass |
| 49 | + |
| 50 | + |
| 51 | +def TestDatasetLPDAAC(): |
| 52 | + dataset = Dataset() |
| 53 | + dataset.getProduct('MOD13Q1') |
| 54 | + dataset.source.download('', ['h12v08']) |
| 55 | + |
| 56 | + |
| 57 | +class LPDAAC(HTTPSource): |
| 58 | + url = 'http://e4ftl01.cr.usgs.gov' |
| 59 | + |
| 60 | + |
| 61 | +class NSIDC(FTPSource): |
| 62 | + url = 'n5eil01u.ecs.nsidc.org' |
| 63 | + host = 'n5eil01u.ecs.nsidc.org' |
| 64 | + user = 'anonymous' |
| 65 | + password = 'anon' |
| 66 | + |
| 67 | + def __init__(self): |
| 68 | + super(NSIDC, self).__init__(self.host) |
| 69 | + self.connect() |
| 70 | + |
| 71 | + |
| 72 | +class MODISHDF(object): |
| 73 | + source = None |
| 74 | + product = dict() #downloadInfo |
| 75 | + |
| 76 | + taskStart = datetime.strptime('1900-01-01', '%Y-%m-%d') |
| 77 | + taskEnd = datetime.strptime('3000-12-31', '%Y-%m-%d') |
| 78 | + geom = None |
| 79 | + geomType = None |
| 80 | + |
| 81 | + directories = [] #dates |
| 82 | + tiles = None |
| 83 | + dataset = None |
| 84 | + bands = None |
| 85 | + |
| 86 | + downloadPath = None |
| 87 | + publishPath = os.getcwd() |
| 88 | + outputs = dict() |
| 89 | + |
| 90 | + def addTask(self, dataset, geom, publishPath, start=None, end=None, qualityValue=None, qualityBand=None, epsg=None, resample=None, source=None, format=None): |
| 91 | + pass |
| 92 | + |
| 93 | + def __init__(self, task, dataset, source): |
| 94 | + LOGGER.debug('LPDAAC init') |
| 95 | + self.dataset = dataset |
| 96 | + self.source = locals()[source]() |
| 97 | + if 'start' in task != "": |
| 98 | + self.taskStart = datetime.strptime(task['start'], '%Y-%m-%d') |
| 99 | + if 'end' in task != "": |
| 100 | + self.taskEnd = datetime.strptime(task['end'], '%Y-%m-%d') |
| 101 | + if not 'geom' in task: |
| 102 | + LOGGER.error('No Geometry given!') |
| 103 | + raise Exception('No Geometry given!') |
| 104 | + sys.exit(1) |
| 105 | + if not 'dataset' in task: |
| 106 | + LOGGER.error('No Dataset given!') |
| 107 | + raise Exception('No Dataset given!') |
| 108 | + sys.exit(1) |
| 109 | + if 'publishPath' not in task: |
| 110 | + task['publishPath'] = os.getcwd() |
| 111 | + self.setPublishPath(task['publishPath']) |
| 112 | + |
| 113 | + self.geom = task['geom'] |
| 114 | + if 'POINT' in self.geom: |
| 115 | + self.geomType = 'POINT' |
| 116 | + else: |
| 117 | + self.geomType = 'POLYGON' |
| 118 | + self.product = task['dataset'] |
| 119 | + self.task = task |
| 120 | + self.bands = dataset.getBands() |
| 121 | + |
| 122 | + self.prepare() |
| 123 | + |
| 124 | + def setPublishPath(self, path): |
| 125 | + self.publishPath = path |
| 126 | + if not os.path.exists(path): os.makedirs(path) |
| 127 | + if not os.path.exists(path+'/data'): os.makedirs(path+'/data') |
| 128 | + |
| 129 | + def setDownloadPath(self, path): |
| 130 | + self.downloadPath = path |
| 131 | + if not os.path.exists(path): os.makedirs(path) |
| 132 | + |
| 133 | + def prepare(self): |
| 134 | + tilesObj = Tiles() |
| 135 | + self.task['tiles'] = tilesObj.getTilesFromGeometry(self.task['geom']) |
| 136 | + self.tiles = '('+'|'.join(self.task['tiles'])+')' |
| 137 | + dates = self.source.listDirectories(self.product['directory']) |
| 138 | + self.directories = [] |
| 139 | + self.dates = [] |
| 140 | + for date in dates: |
| 141 | + dt = datetime.strptime(date, '%Y.%m.%d') |
| 142 | + if dt >= self.taskStart and dt <= self.taskEnd: |
| 143 | + self.directories.append(date) |
| 144 | + self.dates.append(dt) |
| 145 | + |
| 146 | + return dict(tiles=self.task['tiles'], dates=self.dates) |
| 147 | + |
| 148 | + def process(self, files, format, epsg=None, resample=None, compress=False): |
| 149 | + #extract dataset |
| 150 | + LOGGER.info('extract bands') |
| 151 | + raster = processing.MODISHDFProcessor(None, self.bands, self.dataset.rastertype, self.publishPath) |
| 152 | + for file in files: |
| 153 | + raster.extractBands(file) |
| 154 | + |
| 155 | + #merge, clip, reproject, compress based on geometry |
| 156 | + LOGGER.info('merge, clip, reproject, compress') |
| 157 | + if self.geomType == 'POINT': |
| 158 | + outputs = raster.clipPoint(self.geom) |
| 159 | + if 'qualityBand' in self.task and 'qualityValue' in self.task: |
| 160 | + qaValue = outputs[self.task['qualityBand']] |
| 161 | + result = raster.processQualityPoint(qaValue, self.task['qualityBand'], self.task['qualityValue']) |
| 162 | + if result == 0: |
| 163 | + for key, band in self.bands: |
| 164 | + if band['imagetype'] == 'qualityInformation': |
| 165 | + continue |
| 166 | + outputs[key] = self.dataset.nodata |
| 167 | + |
| 168 | + elif self.geomType == 'POLYGON': |
| 169 | + outputs = raster.clipPolygon(self.geom, format, epsg, resample) |
| 170 | + if 'qualityBand' in self.task and 'qualityValue' in self.task: |
| 171 | + qualityOutputs = raster.processQuality(self.task['qualityBand'], self.task['qualityValue']) |
| 172 | + for key in qualityOutputs: |
| 173 | + outputs[key] = qualityOutputs[key] |
| 174 | + else: |
| 175 | + LOGGER.error('No valid geometry found -> no clipping!') |
| 176 | + #raise Exception('No valid geometry found!') |
| 177 | + |
| 178 | + return outputs |
| 179 | + |
| 180 | + def setOutputs(self, files, date): |
| 181 | + if self.geomType == 'POINT': |
| 182 | + for key in files: |
| 183 | + self.outputs[key].addEntry(date, files[key]) |
| 184 | + else: |
| 185 | + if not os.path.exists(self.publishPath+'/data/files'): os.makedirs(self.publishPath+'/data/files') |
| 186 | + for key in files: |
| 187 | + LOGGER.info('move file '+files[key]+' to '+self.publishPath+'/data/files/') |
| 188 | + shutil.move(files[key], self.publishPath+'/data/files/') |
| 189 | + self.outputs[key].addEntry(date, self.publishPath+'/data/files/'+os.path.basename(files[key])) |
| 190 | + |
| 191 | + def download(self, date, path=""): |
| 192 | + if path == "": |
| 193 | + if self.downloadPath != None: |
| 194 | + downloadPath = self.downloadPath+'/'+self.dataset.product+'/'+date |
| 195 | + else: |
| 196 | + downloadPath = self.publishPath+'/raw/'+date |
| 197 | + else: |
| 198 | + downloadPath = path |
| 199 | + if not os.path.exists(downloadPath): os.makedirs(downloadPath) |
| 200 | + |
| 201 | + self.source.listFiles(self.product['directory']+'/'+date) |
| 202 | + self.source.filterFiles('.*\.'+self.tiles+'\..*\.hdf$') |
| 203 | + localFiles = self.source.downloadFiles(downloadPath) |
| 204 | + return localFiles |
| 205 | + |
| 206 | + def mergeOutputs(self, format="VRT", fileformat="HDF4Image"): |
| 207 | + datapath = self.publishPath+'/data/' |
| 208 | + extension = processing.GDAL().getFileExtension(format) |
| 209 | + for key, band in self.bands.items(): |
| 210 | + datafilter = datapath+'files/'+self.dataset.shortname+'.*.'+key+'.*'+processing.GDAL().getFileExtension(fileformat) |
| 211 | + output = datapath+self.dataset.shortname+'.'+key+extension |
| 212 | + if format == 'VRT': |
| 213 | + processing.GDAL().gdalbuildvrt(datafilter, output, True) |
| 214 | + else: |
| 215 | + processing.GDAL().gdal_merge(datafilter, output, format) |
| 216 | + self.outputs[key].setMergedFile(output) |
| 217 | + |
| 218 | + def buildTSRasterOutput(self): |
| 219 | + LOGGER.info('TODO: build time-series raster output') |
| 220 | + pass |
| 221 | + |
| 222 | + def ingest(self, format): |
| 223 | + # set output objects |
| 224 | + if self.geomType == 'POINT': |
| 225 | + from pyEOM.outputs import sensor |
| 226 | + for key, band in self.bands.items(): |
| 227 | + self.outputs[key] = sensor.Sensor(key, self.dataset.shortname, self.publishPath, dict(geom=self.geom, nodata=band['nodata'], scale=band['scale'], offset=band['offset'])) |
| 228 | + else: |
| 229 | + from pyEOM.outputs import raster |
| 230 | + for key, band in self.bands.items(): |
| 231 | + self.outputs[key] = raster.TSRaster(key, self.dataset.shortname, self.publishPath, dict(geom=self.geom, nodata=band['nodata'], scale=band['scale'], offset=band['offset'])) |
| 232 | + |
| 233 | + for date in self.directories: |
| 234 | + files = self.download(date) |
| 235 | + outputs = self.process(files, format) |
| 236 | + self.setOutputs(outputs, date) |
| 237 | + |
| 238 | + if self.geomType == 'POLYGON': |
| 239 | + self.mergeOutputs('VRT', format) |
| 240 | + |
| 241 | + return self.outputs |
| 242 | + |
| 243 | + |
0 commit comments