Array/Image Interpolation in python

Question

I have 365 2d numpy arrays for the every day of the year, displaying an image like this:

I have them all stacked in a 3d numpy array. Pixels with a value that represents cloud i want to get rid of, i want to search through the previous 7 days, or next 7 days (previous 7 layers, next 7 layers) to find a value other than cloud, and then replace the cloud value with other possible values for that pixel (values experienced in the other days/layers for the corresponding pixel).

I am new to python, and am a bit lost.

Any ideas?

Thanks

"replace the cloud value with the other found value" -- there's too much ambiguity in this definition. Please describe the algorithm in more detail. What's the precedence among future/past days and materials? — Brian Cain
– Brian Cain, Commented Nov 28, 2012 at 15:07

Dunes · Accepted Answer · 2012-11-28 18:36:21Z

3

You are essentially trying to write a filter for your array.

First you need to write a function that when given an array of values, with the middle one being the element currently examined, will return some computation of those values. In your case the function will expect to take 1-d array and returns the element nearest to the middle index that is not cloud:

import numpy as np
from scipy.ndimage.filters import generic_filter

_cloud = -1

def findNearestNonCloud(elements):
    middleIndex = len(elements) / 2
    if elements[middleIndex] != _cloud:
        return elements[middleIndex] # middle value is not cloud

    nonCloudIndices, = np.where(elements != _cloud)
    if len(nonCloudIndices) == 0:
        return elements[middleIndex] # all values were cloud

    prevNonCloudIndex = np.where(nonCloudIndices < middleIndex, 
            nonCloudIndices, -1).max()
    nextNonCloudIndex = -np.where(nonCloudIndices > middleIndex, 
            -nonCloudIndices, 1).min()
    # -1 means no non-cloud index

    # pick index closest to middle index    
    if (abs(prevNonCloudIndex - middleIndex) 
            <= abs(nextNonCloudIndex - middleIndex)):
        return elements[prevNonCloudIndex]
    else:
        return elements[nextNonCloudIndex]

Now you need to apply this function to the elements you're interesting. To do this you'll need a mask that denotes which other elements you'll be interested in with regard to a specific element.

from scipy.ndimage.filters import generic_filter

# creates 5 days worth of a 3x3 plot of land
input = np.ones((5, 3, 3)) * _cloud
input[0,:,:] = 10 # set first "image" to all be 10s
input[4,0,0] = 12 # uppper left corner of fourth image is set to 12
print "input data\n", input, "\n"

mask = (5, 1, 1)
# mask represents looking at the present day, 2 days in the future and 2 days in 
# the past for 5 days in total.

print "result\n", generic_filter(input, findNearestNonCloud, size=mask)
# second and third images should mirror first image, 
# except upper left corner of third image should be 12

answered Nov 28, 2012 at 18:36

Dunes

42.1k7 gold badges87 silver badges107 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

David Colpey Over a year ago

Thanks, will give this a go in the morning. Will get back to you. Thanks again.

Dunes Over a year ago

I put my dimensions as [z,y,x]. Yours may be different.

David Colpey · Accepted Answer · 2012-12-31 12:35:48Z

1

I solved it by this:

interpdata = []
j = 0
for i in stack:
    try:
        temp = np.where( stack[j] == 50, stack[j-1], modis[j] )
        temp = np.where( temp == 50, stack[j+1], temp )
        temp = np.where( temp == 50, stack[j-2], temp )
        temp = np.where( temp == 50, stack[j+2], temp )
        temp = np.where( temp == 50, stack[j-3], temp )
        temp = np.where( temp == 50, stack[j+3], temp ) 
        temp = np.where( temp == 50, stack[j-4], temp )
        temp = np.where( temp == 50, stack[j+4], temp )
    except IndexError:
        print 'IndexError Passed'       
        pass
    else:
        pass
    interpdata [j, :, :] = temp
    j = j + 1

answered Dec 31, 2012 at 12:35

David Colpey

812 silver badges9 bronze badges

Comments

mgilson · Accepted Answer · 2012-11-28 15:28:15Z

0

I would think that you could do something like:

data = somehow_get_your_3d_data() #indexed as [day_of_year,y,x]
for i,dat in enumerate(data):
    weeks2 = data[max(i-7,i):min(i+7,len(data)), ... ]
    new_value = get_new_value(weeks2) #get value from weeks2 here somehow
    dat[dat == cloud_value] = new_value

answered Nov 28, 2012 at 15:28

mgilson

312k70 gold badges656 silver badges722 bronze badges

2 Comments

Avaris Over a year ago

I guess it should be max(i-7, 0).

David Colpey Over a year ago

Okay, But how am i meant to extract the value form the relevant pixel, i.e. 'get_new_value'?

Collectives™ on Stack Overflow

Array/Image Interpolation in python

3 Answers 3

2 Comments

Comments

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

3 Answers 3

2 Comments

Comments

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related