for e in NDVI_crops:

NDVI_crops2 = e.data[~e.mask] # Exclude non-cropland pixel values

NDVI_crops2 = NDVI_crops2[~np.isnan(NDVI_crops2)] # Exclude poor quality pixel values

n = len(NDVI_crops2) # Count of array

min_val = float(format((np.min(NDVI_crops2)), '.4f')) # Minimum value in array

max_val = float(format((np.max(NDVI_crops2)), '.4f')) # Maximum value in array

range_val = (min_val, max_val) # Range of values in array

mean = float(format((np.mean(NDVI_crops2)), '.4f')) # Mean of values in array

std = float(format((np.std(NDVI_crops2)), '.4f')) # Standard deviation of values in array

var = float(format((np.var(NDVI_crops2)), '.4f')) # Variance of values in array

median = float(format((np.median(NDVI_crops2)), '.4f')) # Median of values in array

quartiles = np.percentile(NDVI_crops2, [25, 75]) # 1st (25) & 3rd (75) quartiles of values in array

upper_quartile = float(format((quartiles[1]), '.4f'))

lower_quartile = float(format((quartiles[0]), '.4f'))

iqr = quartiles[1] - quartiles[0] # Interquartile range

iqr_upper = upper_quartile + 1.5 * iqr # 1.5 IQR of the upper quartile

iqr_lower = lower_quartile - 1.5 * iqr # 1.5 IQR of the lower quartile

top = float(format(np.max(NDVI_crops2[NDVI_crops2 =iqr_lower]),'.4f')) # Lowest datum within 1.5 IQR of lower quartile