From 09f142644e45e3d95d6b1910d9fdb7dc58a1790e Mon Sep 17 00:00:00 2001 From: "Bobby H. Braswell" Date: Tue, 19 Nov 2019 21:58:57 -0500 Subject: [PATCH] new new features because south america --- multitemporal/bin/shifttime.pyx | 45 +++++++++++++++++++++++++++++++++ multitemporal/bin/trimyr.pyx | 41 ++++++++++++++++++++++++++++++ multitemporal/mt.py | 39 ++++++++++++++++++++++------ multitemporal/test/test_1.py | 27 ++++++++++++++++++++ 4 files changed, 144 insertions(+), 8 deletions(-) create mode 100644 multitemporal/bin/shifttime.pyx create mode 100644 multitemporal/bin/trimyr.pyx diff --git a/multitemporal/bin/shifttime.pyx b/multitemporal/bin/shifttime.pyx new file mode 100644 index 0000000..612d605 --- /dev/null +++ b/multitemporal/bin/shifttime.pyx @@ -0,0 +1,45 @@ +import numpy as np +cimport numpy as np +cimport cython + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) + +def get_nout(int nin, np.ndarray[np.float32_t, ndim=1, negative_indices=False] params not None): + return nin + +def get_nyrout(int nyr, np.ndarray[np.float32_t, ndim=1, negative_indices=False] params not None): + return nyr + +def shifttime(np.ndarray[np.float32_t, ndim=3, negative_indices=False] data not None, + float missingval, + np.ndarray[np.float32_t, ndim=1, negative_indices=False] params not None): + + cdef unsigned int nfr = data.shape[0] + cdef unsigned int nyr = data.shape[1] + cdef unsigned int npx = data.shape[2] + + cdef int offset = params[0] - 1 + + cdef np.ndarray[np.float32_t, ndim=3] result = np.zeros((nfr, nyr, npx), dtype='float32') + + cdef int i, j, k + cdef int t, nt + cdef int t1, i1, j1 + cdef float count, ave + + nt = nfr*nyr + + for k in range(npx): + for t in range(nt): + t1 = t - offset + if t1 < 0: + continue + i = t % nfr + j = t / nfr + i1 = t1 % nfr + j1 = t1 / nfr + result[i1,j1,k] = data[i,j,k] + + return result diff --git a/multitemporal/bin/trimyr.pyx b/multitemporal/bin/trimyr.pyx new file mode 100644 index 0000000..d9e5be2 --- /dev/null +++ b/multitemporal/bin/trimyr.pyx @@ -0,0 +1,41 @@ +import numpy as np +cimport numpy as np +cimport cython + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) + +def get_nout(int nin, np.ndarray[np.float32_t, ndim=1, negative_indices=False] params not None): + return nin + +def get_nyrout(int nyr, np.ndarray[np.float32_t, ndim=1, negative_indices=False] params not None): + return params[0] - params[1] + 1 + +def trimyr(np.ndarray[np.float32_t, ndim=3, negative_indices=False] data not None, + float missingval, + np.ndarray[np.float32_t, ndim=1, negative_indices=False] params not None): + + cdef unsigned int nfr = data.shape[0] + cdef unsigned int nyr = data.shape[1] + cdef unsigned int npx = data.shape[2] + + cdef int nout = get_nout(nfr, params) + cdef int nyrout = get_nyrout(nyr, params) + + cdef int yr1 = params[0] - 1 + cdef int yr2 = params[1] - 1 + + cdef np.ndarray[np.float32_t, ndim=3] result = np.zeros((nout,nyrout,npx), dtype='float32') + + cdef unsigned int i,j,k,idx + + for k in range(npx): + jdx = 0 + for j in range(nyr): + if j>= yr1 and j <= yr2: + for i in range(nfr): + result[i,jdx,k] = data[i,j,k] + jdx = jdx + 1 + + return result diff --git a/multitemporal/mt.py b/multitemporal/mt.py index e33c644..6eb56e7 100755 --- a/multitemporal/mt.py +++ b/multitemporal/mt.py @@ -113,7 +113,7 @@ def worker(shared, job): data[ib, ifr, iyr, wgood] = \ source['offset'] + source['scale']*values[wgood] del fp - + sourcenames = [source['name'] for source in sources] results = {} @@ -131,7 +131,7 @@ def worker(shared, job): d = np.array([results[si] for si in step['inputs']]) if d.shape[0] == 1: - d = d.reshape(d.shape[1], nyr, npx) + d = d.reshape(d.shape[1], d.shape[2], npx) results[step['name']] = step['function'](d, missing_out, step['params']) if step.get('output', False): @@ -182,7 +182,15 @@ def run(projdir, outdir, projname, sources, steps, for i, source_path in enumerate(source_paths): source_bn = os.path.basename(source_path) - datestr = re.findall(source['regexp'], source_bn)[0] + match = re.match(source['regexp'], source_bn) + if len(match.groups()) == 1: + # the normal case of only a variable date string + datestr = match.groups()[0] + else: + # unusual case where you have different locations piled together + datestr = match.group('datestr') + if kwargs.get('tileid', None) != match.group('tileid'): + continue if not ymd: # default: YYYYDDD date = datetime.datetime.strptime(datestr, '%Y%j') @@ -306,7 +314,8 @@ def run(projdir, outdir, projname, sources, steps, for thisinput in step['inputs']: # if this input is in the sources then we know it's a starting point for the pipeline: if thisinput in [source['name'] for source in sources]: - thisnin = nfr # nfr == len(doys) + thisnin = nfr + thisnyrin = nyr step['initial'] = True else: # handle intermediary steps in the pipeline; only one parent is allowed @@ -316,17 +325,26 @@ def run(projdir, outdir, projname, sources, steps, 'specified input: {}\n\tpossible inputs: {}' .format(thisinput, [s['name'] for s in steps])) thisnin = parentsteps[0]['nout'] + thisnyrin = parentsteps[0]['nyrout'] if 'nin' in step: assert step['nin'] == thisnin, "Number of inputs do not match" else: step['nin'] = thisnin + if 'nyrin' in step: + assert step['nyrin'] == thisnyrin, "Number of years do not match" + else: + step['nyrin'] = thisnyrin # set the number of outputs for each step - step['nout'] = int(mod.get_nout(step['nin'], step['params'])) + # TODO: any unexpected behavior here? try: - step['nyrout'] = int(mod.get_nyrout(nyr, step['params'])) + step['nout'] = int(mod.get_nout(step['nin'], step['params'])) except: - step['nyrout'] = nyr + step['nout'] = step['nin'] + try: + step['nyrout'] = int(mod.get_nyrout(step['nyrin'], step['params'])) + except: + step['nyrout'] = step['nyrin'] if step.get('output', False): print("output", mod, (step['nout'], step['nyrout'], height*width)) OUTPUT[step['name']] = sharedmem.empty( @@ -365,10 +383,15 @@ def run(projdir, outdir, projname, sources, steps, prog += 1 print('mt {:0.02f} complete.\r'.format(pct)) results.append(r) - else: + elif nproc == 1: results = [] for job in jobs: results.append(func(job)) + elif nproc < 0: + # secret way to test + results = [func(jobs[-nproc-1])] + else: + raise Exception("nproc can't be zero") # write outputs if not os.path.exists(outdir): diff --git a/multitemporal/test/test_1.py b/multitemporal/test/test_1.py index d75e759..f468da0 100644 --- a/multitemporal/test/test_1.py +++ b/multitemporal/test/test_1.py @@ -8,6 +8,7 @@ from multitemporal import mt + def data_dir(): return os.path.join(os.path.dirname(__file__), 'data') @@ -36,6 +37,32 @@ def test_passthrough(nproc, tmpdir): # both pytest & numpy have approximate-equality abilities if needed assert (expected == actual).all() + +test_passthrough2_args = { + 'compthresh': 0.01, # so smaller dataset will work + "dperframe": 1, + "sources" : [{"name": "ndvi", "regexp": "^(?P\\d{7})_L.._ndvi-(?P\w{3}).tif$", "bandnum": 1}], + "steps" : [{"module": "passthrough", "params": [], "inputs": ["ndvi"], "output": True}]} + + +@pytest.mark.parametrize("nproc", [1, 2]) +def test_passthrough2(nproc, tmpdir): + """Use the passthrough module as a way to test mt throughput, this time with weird regex""" + # refactor out for additional tests and follow pattern in files: + input_dir = os.path.join(data_dir(), 'input') + output_bn = 'tpt_proj_passthrough.tif' + expected_fp = os.path.join(data_dir(), 'expected', output_bn) + actual_fp = str(tmpdir.join(output_bn)) + + mt.run(projname='tpt_proj', projdir=input_dir, outdir=str(tmpdir), nproc=nproc, tileid='toa', + **test_passthrough2_args) + + actual = gdal.Open(actual_fp).ReadAsArray() + expected = gdal.Open(expected_fp).ReadAsArray() + # both pytest & numpy have approximate-equality abilities if needed + assert (expected == actual).all() + + def test_two_sources(tmpdir): """As test_passthrough, but with two sources.