2

I'm working with a time series stack of 28 raster bands. Each band is an image for a specific date, spanning all of the year 2020, and they are stacked in order. The stack (change_point_raster) is in the form of a 3D array where the first dimension is the date band. I also have an associated raster (change_point_index) that is just one band, where each pixel represents the number of the band in change_point_raster where there is detected change. I want to do some statistics before and after the change occurred for each pixel. To do this, I want to get the bands before and after the change band, then find the ratio between the two band values for each pixel in the stack. I'm doing this in Python. I'm kind of a beginner in python, so I'm unsure what to do, but here is a code snippet of what I have:

#Identify change date (or band number) and subset/index time series stack by before change and after change.

cpr_before = next(x for x, val in enumerate(change_point_raster)
                                  if val < change_point_index)

cpr_after = next(x for x, val in enumerate(change_point_raster)
                                  if val > change_point_index)

I'm getting an error that says "ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()". It seems like the problem is that I have a raster stack instead of just one raster, but I'm unsure how to solve this with a loop, or any other method.

3
  • Don't you need square brackets around list comprehensions? – wingnut May 2 at 0:44
  • @wingnut yes, but... this is not a list comprehension, it's a generator stackoverflow.com/q/364802/737471 – user2856 May 2 at 7:04
  • Oh yeah, I see it now. It returns a 2D array. – wingnut May 2 at 7:12
0

Try:

cpr_before = next(x for val, x in enumerate(change_point_raster)
                              if val < change_point_index)

Instead of

cpr_before = next(x for x, val in enumerate(change_point_raster)
                              if val < change_point_index)

Ithink you want to return an array based on an index, not the other way around. Same for cpr_after.

The following code creates dummy data with dummy index breaks. Then it looks up the index table for each pixel, and calculates the mean on 'days' before the break, and on 'days' after the break. It's simplified but should help.

# imports
import numpy as np
import matplotlib.pyplot as plt

# simulate data
w = 30
h = 30
N = 10

# random indices
index = np.random.randint(0,N,[w,h])

# make the data look different before and after the index
data = np.zeros([N,w,h])
for i in range(w):
    for j in range(h):
        data[:index[i,j],:,:] = np.random.randn(index[i,j],w,h)
        data[index[i,j]:,:,:] = np.random.randn(N-index[i,j],w,h)+3

# Analysis
# at each pixel, get the data before and after the break specified by the index layer
# and save the means before and after to separate layers
arr_before = np.zeros([w,h])
arr_after = np.zeros([w,h])
for i in range(w):
    for j in range(h):
        if index[i,j]==0:
            arr_before[i,j] = data[0,i,j]
        else:
            arr_before[i,j] = np.mean(data[:index[i,j],i,j])
        if index[i,j] == N-1:
            arr_after[i,j] = data[-1,i,j]
        else:
            arr_after[i,j]  = np.mean(data[index[i,j]:,i,j])


# Plot means before and after the index break
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.imshow(arr_before, vmin=-5, vmax=5)
plt.colorbar()
plt.title('Before')
plt.subplot(1,2,2)
plt.imshow(arr_after, vmin=-5, vmax=5)
plt.colorbar()
plt.title('After')
3
  • thanks, I think you're right. However, when I try that I'm still getting this error: ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all() – Katie Walker May 2 at 22:09
  • It may be because change_point_raster is a 3D array and change_point_index is a 2D array? – Katie Walker May 2 at 22:16
  • Yes, I added some example code in the answer to do it pixel by pixel. It calculates means of synthetic data before and after random index breaks and then plots their rasters. You'll need to modify it for your needs, but it should be easy. – wingnut May 3 at 5:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.