I am trying to mask an array and then add the array to a list, so I can
then go on and calculate the max, min and mean of that list.
The mask seems to work when I check each array. I check each array by
finding the max, mean and mean and comparing with the unmasked array (they
are different).
However, when I add the array to the list and then check the max, min and
mean and compare to the list of unmasked array I find that they are the
same (which they shouldn't be).
The problem seems to be occurring when I add the masked array to the list.
Is there something I am doing wrong or not doing??
Any feedback will be greatly appreciated.
import numpy as np
import matplotlib.pyplot as plt
from numpy import ma as MA
from mpl_toolkits.basemap import Basemap
from datetime import datetime
import os
from StringIO import StringIO
from osgeo import gdal, gdalnumeric, ogr, osr
import glob
from datetime import date, timedelta
import matplotlib.dates as mdates
import time
shapefile=r"d:/Vic/Vic_dissolve.shp"
## Create masked array from shapefile
xmin,ymin,xmax,ymax=[111.975,-9.975, 156.275,-44.525]
ncols,nrows=[886, 691]
maskvalue = 1
xres=(xmax-xmin)/float(ncols)
yres=(ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)
0
src_ds = ogr.Open(shapefile)
src_lyr=src_ds.GetLayer()
dst_ds = gdal.GetDriverByName('MEM').Create('',ncols, nrows, 1
,gdal.GDT_Byte)
dst_rb = dst_ds.GetRasterBand(1)
dst_rb.Fill(0) #initialise raster with zeros
dst_rb.SetNoDataValue(0)
dst_ds.SetGeoTransform(geotransform)
err = gdal.RasterizeLayer(dst_ds, [maskvalue], src_lyr)
dst_ds.FlushCache()
mask_arr=dst_ds.GetRasterBand(1).ReadAsArray()
np.set_printoptions(threshold='nan')
mask_arr[mask_arr == 255] = 1
newmask=MA.masked_equal(mask_arr,0)
### calculate monthly summary stats for VIC Only
rainmax=[]
rainmin=[]
rainmean=[]
rainmaxaust=[]
rainminaust=[]
rainmeanaust=[]
yearmonthlist=[]
yearmonth_int=[]
GLOBTEMPLATE = r"e:/Rainfall/rainfall-{year}/r{year}{month:02}??.txt"
def accumulate_month(year, month):
files = glob.glob(GLOBTEMPLATE.format(year=year, month=month))
monthlyrain=[]
monthlyrainaust=[]
for ifile in files:
f=np.genfromtxt(ifile,skip_header=6)
viconly_f=np.ma.masked_array(f, mask=newmask.mask)
#print "f:", f.max(), f.mean()
#print "viconly_f:", viconly_f.max(), viconly_f.mean()
monthlyrain.append(viconly_f)
monthlyrainaust.append(f)
yearmonth=str(year)+str(month)
yearmonthlist.append(yearmonth)
r_max, r_mean, r_min=np.max(monthlyrain), np.mean(monthlyrain),
np.min(monthlyrain)
ra_max, ra_mean, ra_min=np.max(monthlyrainaust),
np.mean(monthlyrainaust), np.min(monthlyrainaust)
rainmax.append(r_max)
rainmean.append(r_mean)
rainmean.append(r_min)
rainmaxaust.append(ra_max)
rainminaust.append(ra_min)
rainmeanaust.append(ra_mean)
print "viconly:", yearmonth,r_max, r_mean, r_min
print "aust:", yearmonth,ra_max, ra_mean, ra_min
###loop through months and years
stop_month = datetime(2011, 10, 01)
month = datetime(2011, 01, 01)
while month < stop_month:
accumulate_month(month.year, month.month)
month += timedelta(days=32)
month = month.replace(day=01)
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
