def __init__(self, filesearch):
     """docstring for init"""
     self.files=glob.glob(filesearch)
     self.files.sort()
     
     if len(self.files)==0:
         raise ValueError("\nERROR: no files match search string:"+filesearch+"\n")
     if verbose: print("Reading:{}".format(self.files[self.curfile]))
     self.tdata  = mygis.read_nc(self.files[self.curfile], self.t_var).data-273.15
     self.rrdata = mygis.read_nc(self.files[self.curfile],self.rr_var).data
     self.acdata = mygis.read_nc(self.files[self.curfile],self.ac_var).data
     self.winddata = mygis.read_nc(self.files[self.curfile],self.u_var).data**2
     self.winddata+= mygis.read_nc(self.files[self.curfile],self.v_var).data**2
     self.winddata=np.sqrt(self.winddata)
     self.husdata = mygis.read_nc(self.files[self.curfile],self.hus_var).data
     self.swdata = mygis.read_nc(self.files[self.curfile],self.sw_var).data
     self.lwdata = mygis.read_nc(self.files[self.curfile],self.lw_var).data
     self.acsnowdata = mygis.read_nc(self.files[self.curfile],self.snow_var).data
     self.snowdata=np.zeros(self.acsnowdata.shape)
     self.snowdata[0]=self.acsnowdata[0]
     self.snowdata[1:]=np.diff(self.acsnowdata,axis=0)
     for i in range(self.snowdata.shape[0]):
         if self.snowdata[i].min()<0:
             self.snowdata[i]=self.acsnowdata[i]
     self.lastsnow=self.acsnowdata[-1]
     
     self.mask = mygis.read_nc(mask_file,"LANDMASK").data[0]==0

def load_wrf_data(stat,time):
    """Load WRF data and calculate the relevant statistic"""
    output_wrf_file=stat+"_"+time+".nc"
    try:
        current=mygis.read_nc(wrfloc+"current_"+output_wrf_file).data
        future=mygis.read_nc(wrfloc+"future_"+output_wrf_file).data
        print("WRF data loaded from pre-calculated summary files")
    except:
        print("    Reading raw data...")
        if time=="annual":
            current=np.concatenate(mygis.read_files(wrfloc+"daily_NARR*"))
            future=np.concatenate(mygis.read_files(wrfloc+"daily_PGW*"))
        else:
            month=time[-2:]
            current=np.concatenate(mygis.read_files(wrfloc+"daily_NARR*"+month+".nc"))
            future=np.concatenate(mygis.read_files(wrfloc+"daily_PGW*"+month+".nc"))
        print("    Generating statistics")
        if stat=="MAP":
            ndays=calc_ndays(time)
            current=current.mean(axis=0)*ndays
            future=future.mean(axis=0)*ndays
        
        elif stat=="wetfrac":
            current=calc_wetfrac(current)
            future=calc_wetfrac(future)
        else:
            raise KeyError("stat not created for WRF:"+stat)
        print("    Writing stat for future reference")
        mygis.write(wrfloc+"current_"+output_wrf_file,current)
        mygis.write(wrfloc+"future_"+output_wrf_file,future)
        
    # if stat=="MAP":
    #     current[current<1e-4]=1e-4
    return (future-current),current

def main (filename, outputfile):
    
    icar_data=ICAR_Reader(filename)
    if verbose:print(icar_data.files[0],icar_data.files[-1])
    
    raindata=[]
    tmindata=[]
    tmaxdata=[]
    tavedata=[]
    if verbose:print("Looping through data")
    for data in icar_data:
        raindata.append(data[0])
        tmindata.append(data[1])
        tmaxdata.append(data[2])
        tavedata.append(data[3])
    
    latvar.data=mygis.read_nc(icar_data.files[0],"lat").data
    lonvar.data=mygis.read_nc(icar_data.files[0],"lon").data
    
    dates_are_mjd = (icar_data.files[0].split("/")[1] == "erai")
    if verbose:print("using modified julian day dates="+str(dates_are_mjd))
    year = int(icar_data.files[0].split("/")[-1].split("_")[1])
    if dates_are_mjd:
        timevar.attributes = era_time_atts
        start_date = date_fun.date2mjd(year,01,01,00,00)
    else:
        start_date = (year-1800) * 365
    timevar.data = np.arange(start_date,start_date+len(raindata))
    if verbose:print(year, timevar.data[0], len(raindata))
    
    if verbose:print("Writing data")
    # write_file(outputfile+"_rain.nc",raindata, varname="precipitation_amount",      varatts=prec_atts)
    # write_file(outputfile+"_tmin.nc",tmindata, varname="daily_minimum_temperature", varatts=tmin_atts)
    # write_file(outputfile+"_tmax.nc",tmaxdata, varname="daily_maximum_temperature", varatts=tmax_atts)
    write_file(outputfile+"_tave",tavedata, varname="daily_average_temperature", varatts=tave_atts)

def load_erai_means(wind_option):
    """docstring for load_erai_means"""
    eraid="erai/"
    varlist=["p","rh","ta","ua","va","z"]
    if (wind_option=="nowind"):
        varlist=["p","rh","ta","z"]
        
    outputdata=[]
    month_mid_point_doy=(start_day_per_month[1:]+np.array(start_day_per_month[:-1]))*0.5
    
    for month in range(1,13):
        curoutput=Bunch(doy=month_mid_point_doy[month-1])
        for v in varlist:
            if v=="ta":
                ta=mygis.read_nc(eraid+"regridded_ERAi_to_cesm_month{0:02}.nc".format(month),v).data
            else:
                curoutput[v]=mygis.read_nc(eraid+"regridded_ERAi_to_cesm_month{0:02}.nc".format(month),v).data
        curoutput["theta"] = ta / units.exner(curoutput.p)
        
        # erai is probably "upside down" so reverse the vertical dimension of all variables
        if curoutput.p[0,0,0]<curoutput.p[-1,0,0]:
            for v in curoutput.keys():
                if v!="doy":
                    curoutput[v]=curoutput[v][::-1,:,:]
        outputdata.append(curoutput)
    return outputdata

    def nextNN(self):
        curfile=self._curfile
        if curfile>=len(self._filenames):
            raise StopIteration
        minx=self.x.min()
        maxx=self.x.max()+1
        miny=self.y.min()
        maxy=self.y.max()+1
        curx=self.x-minx
        cury=self.y-miny
        data=list()
        for thisvar in self._vars:
            ncfile=mygis.read_nc(self._filenames[curfile],var=thisvar,returnNCvar=True)
            if self.npos==1:
                data.append(ncfile.data[miny:maxy,minx:maxx][cury,curx])
            else:
                data.append(ncfile.data[self.posinfile,miny:maxy,minx:maxx][cury,curx])
            ncfile.ncfile.close()
        self.posinfile+=1

        if (self.npos==1) or (self.posinfile>=self.npos):
            self._curfile+=1
            self.posinfile=0
            if self.npos>1:
                if self._curfile<len(self._filenames):
                    d=mygis.read_nc(self._filenames[self._curfile],var=self._vars[0],returnNCvar=True)
                    ntimes=d.data.shape[0]
                    d.ncfile.close()
                    self.npos=ntimes
            
        
        return data

def main(start_date="19900101"):
    """convert a file (or files) from downscaling code to maps"""
    files=find_files(start_date)
    files.sort()
    geo=load_geo(global_basefile)
    times=mygis.read_nc(files[0],"time").data
    ntimes=times.size
    
    tmp=mygis.read_nc(files[0],"coefficient",returnNCvar=True)
    output_data=np.zeros((tmp.data.shape[1],ntimes,geo.lon.shape[0],geo.lon.shape[1]))
    tmp.ncfile.close()
    print(output_data.shape)
    for f in files:
        print(f)
        data=mygis.read_nc(f,"coefficient").data
        locations=get_xy(f,geo)
        for i in range(len(locations.x)):
            output_data[:,:,locations.y[i],locations.x[i]]=data[0,:,:,i]
        
    print("Writing output file")
    mygis.write(global_output_file,output_data,varname="coefficient",dtype="d",dims=("variable","time","latitude","longitude"),
                    extravars=[ Bunch(data=times,name="time",dims=("time",),dtype="d",
                                      attributes=Bunch(units="seconds since 1970-01-01 00:00:00.0 0:00")),
                                Bunch(data=geo.lat,name="latitude",dims=("latitude","longitude"),dtype="f",
                                      attributes=Bunch(units="degrees")),
                                Bunch(data=geo.lon,name="longitude",dims=("latitude","longitude"),dtype="f",
                                      attributes=Bunch(units="degrees"))])

def load_base_data(swefile="SWE_daily.nc",info="4km_wrf_output.nc",res=4,year=5):
    # wrf.load_base_data(swefile="wrfout_d01_2008-05-01_00:00:00",res=2)
    
    if res==2:
        forest=[11,12,13,14,18]
        exposed=[1,2,3,4,5,7,8,9,10,16,17,19,20,22,23,24,25,26,27,15,21] #warning 15 and 21 sound like they should be forest, but on the map they are in open areas
        mayday=0
        info=swefile
        dz=100
    else:
        forest=[1,5]
        exposed=[7,10]
        mayday=212
        for i in range(year):
            mayday+=365
        mayday+=np.floor(year/4)
        
    print("Loading data")
    vegclass=myio.read_nc(info,"IVGTYP").data[0,...]
    mask=np.zeros(vegclass.shape)
    for f in forest:
        forested=np.where(vegclass==f)
        mask[forested]=1
    for e in exposed:
        exposed=np.where(vegclass==e)
        mask[exposed]=2
    lat=myio.read_nc(info,"XLAT").data[0,...]
    lon=myio.read_nc(info,"XLONG").data[0,...]
    dem=myio.read_nc(info,"HGT").data[0,...]
    all_snow=myio.read_nc(swefile,"SNOW").data/1000.0
    snow=all_snow[mayday,:,:]
    data_by_years=[all_snow[baseday:baseday+365,:,:] for baseday in range(0,all_snow.shape[0]-20,365)]
    return Bunch(data=snow, lat=lat,lon=lon, dem=dem, lc=mask, yearly=data_by_years)

def load_geo(filename):
    """load a 2d lat and lon grid from filename"""
    latnames=["lat","latitude","XLAT"]
    lonnames=["lon","longitude","XLONG"]
    lat=None
    lon=None
    
    for l in latnames:
        if lat==None:
            try:
                lat=mygis.read_nc(filename,l).data
            except:
                lat=None
    for l in lonnames:
        if lon==None:
            try:
                lon=mygis.read_nc(filename,l).data
            except:
                lon=None

    ymin=np.where(lat>global_bounds.lat.min)[0][0]
    ymax=np.where(lat>global_bounds.lat.max)[0][0]
    lat=lat[ymin:ymax]
    
    xmin=np.where(lon>global_bounds.lon.min)[0][0]
    xmax=np.where(lon>global_bounds.lon.max)[0][0]
    lon=lon[xmin:xmax]

    if len(lon.shape)==1:
        lon,lat=np.meshgrid(lon,lat)
        
    return Bunch(lon=lon,lat=lat)

def load_elev_comparison(swefile="SWE_Daily0600UTC_WesternUS_2010.dat",outputfile="snodas_by_elev.png",domainname="Front Range +",domain=None,dz=100):
    import wsc.compare2lidar as c2l

    demfile="snodas_dem.nc"
    forestfile="snodas_forest.nc"
    forest=[1]
    bare=[0]

    mayday=120
    if domain is None:
        minx=1800; miny=600; maxx=None;maxy=1100
    else:
        miny=domain[0]; maxy=domain[1]; minx=domain[2]; maxx=domain[3]

    print("Loading data")
    vegclass=myio.read_nc(forestfile).data[miny:maxy,minx:maxx]
    forested=np.where(vegclass==forest[0])
    exposed=np.where(vegclass==bare[0])
    mask=np.zeros(vegclass.shape)
    mask[forested]=1
    mask[exposed]=2

    dem=myio.read_nc(demfile).data[miny:maxy,minx:maxx]
    snodas=load(swefile)
    snow=snodas.data[mayday,miny:maxy,minx:maxx]

    print("Binning")
    banded=c2l.bin_by_elevation(snow,dem,mask,dz=dz)

    print("Plotting")
    c2l.plot_elevation_bands(banded,outputfile,title="SNODAS SWE over "+domainname)

def write_monthly_pgw(current,future,geofile):
    """write out a file containing the ratios future/current for each month and lat/lon data"""
    
    print("writing PGW results")
    if esgvarname=="pr":
        output_data=future/current
        outputfilename="PGW_ratio_file.nc"
        data_atts=Bunch(long_name="Ratio between future:current monthly precipitaton", units="mm/mm")
    else:
        output_data=future-current
        outputfilename="PGW_difference_file.nc"
        data_atts=Bunch(long_name="Difference between future-current monthly Temperature", units="K")
        
    lat=io.read_nc(geofile,"lat").data
    lon=io.read_nc(geofile,"lon").data
    time=np.arange(12)
    
    lat_atts =Bunch(long_name="latitude",  units="degrees")
    lon_atts =Bunch(long_name="longitude", units="degrees")
    time_atts=Bunch(long_name="Month of the year",units="month",description="0=January,11=December,etc.")
    
    datadims=("time","lat","lon")
    evars=[ Bunch(data=lat, name="lat", dtype="f",attributes=lat_atts,dims=("lat",)),
            Bunch(data=lon, name="lon", dtype="f",attributes=lon_atts,dims=("lon",)),
            Bunch(data=time,name="time",dtype="f",attributes=time_atts,dims=("time",)),
            ]
    
    io.write(outputfilename,output_data,dtype="f",varname="data",dims=datadims,attributes=data_atts,extravars=evars)

def load_tskin(filename,tsvarname, maskvarname):
    """load skin temp data from wrf output file: filename
    
    tskin uses it's own procedure so that we can average over the past 10 days for land points
    ideally it will create a higher resolution output based on a high-resolution land-sea mask
    """
    if verbose:print("Loading :"+maskvarname)
    ncmask=mygis.read_nc(filename,maskvarname,returnNCvar=True)
    mask=ncmask.data[0]
    ncmask.ncfile.close()
    land=np.where(mask==1)
    
    if verbose:print("Loading :"+tsvarname)
    
    ncts=mygis.read_nc(filename,tsvarname,returnNCvar=True)
    atts=mygis.read_atts(filename,tsvarname)
    dims=ncts.data.dimensions
    
    if verbose:print("    "+str(dims))
    if verbose:print("    "+str(atts))
    
    data=ncts.data[:]
    ncts.ncfile.close()
    
    for i in range(data.shape[0]-1,0,-1):
        start=max(i-(10*steps_per_day),0)
        data[i][land]=data[start:i+1].mean(axis=0)[land]
    
    return Bunch(name=tsvarname, data=data, attributes=atts, dims=dims, dtype='f')

def load_monthly_ar(filesearch,geolut):
    """load monthly means for a given filesearch"""
    files=[]
    for search in filesearch:
        files.extend(glob.glob(search))
    files.sort()
    ncdata=io.read_nc(files[0],esgvarname,returnNCvar=True)
    outputdata=np.zeros((12,ncdata.data.shape[1],ncdata.data.shape[2]))
    ncdata.ncfile.close()
    
    for f in files:
        data=io.read_nc(f,esgvarname).data
        startpoint=0
        for i in range(12):
            endpoint=startpoint+month_lengths[i]
            outputdata[i,...]+=data[startpoint:endpoint,:,:].mean(axis=0)
            startpoint=endpoint
    outputdata/=len(files)
    
    data=np.zeros((12,geolut.shape[0],geolut.shape[1]))
    for i in range(4):
        y=geolut[:,:,i,0].astype('i')
        x=geolut[:,:,i,1].astype('i')
        data+=np.float32(outputdata[:,y,x]*geolut[np.newaxis,:,:,i,2])
    
    return data

def load_hi_res_dem():
    """docstring for load_hi_res_dem"""
    demf=mygis.read_nc(dem_file,"elev_m",returnNCvar=True)
    latf=mygis.read_nc(dem_file,"lat",returnNCvar=True)
    lonf=mygis.read_nc(dem_file,"lon",returnNCvar=True)
    if DEBUG:
        x0=2500;x1=4000;y0=12000;y1=10500
    # x0=0;x1=None;y0=-1;y1=0
        dem=demf.data[y0:y1:-1,x0:x1]
        lat=latf.data[y0:y1:-1]
        lon=lonf.data[x0:x1]
    else:
        dem=demf.data[::-1,:]
        lat=latf.data[::-1]
        lon=lonf.data[:]
    demf.ncfile.close()
    latf.ncfile.close()
    lonf.ncfile.close()
    
    dx=lon[1]-lon[0]
    dy=lat[1]-lat[0]
    nx=len(lon)
    ny=len(lat)
    lon,lat=np.meshgrid(lon,lat)
    
    return Bunch(data=dem,lat=lat,lon=lon,
                 startx=lon[0,0]-dx/2,starty=lat[0,0]-dy/2,
                 dx=dx,dy=dy,nx=nx,ny=ny)

def load_cesm_means(wind_option):
    """docstring for load_cesm_means"""
    cesmd="means/"
    varlist=["p","rh","theta","u","v"]
    if wind_option=="nowind":
        varlist=["p","rh","theta"]
        
    outputdata=[]
    month_mid_point_doy=(start_day_per_month[1:]+np.array(start_day_per_month[:-1]))*0.5
    
    for month in range(1,13):
        curoutput=Bunch(doy=month_mid_point_doy[month-1])
        try:
            for v in varlist:
                curoutput[v]=mygis.read_nc(cesmd+"month{0:02}_mean_{1}.nc".format(month,v),v).data
            curoutput["z"]=mygis.read_nc(cesmd+"annual_mean_z.nc","z").data
        except:
            for v in varlist:
                curoutput[v]=mygis.read_nc("month{0:02}_mean_{1}.nc".format(month,v),v).data
            curoutput["z"]=mygis.read_nc("annual_mean_z.nc","z").data
            
        if curoutput["p"][0,0,0]<1300:
           curoutput["p"][0,0,0]*=100.0 # convert hPa to Pa 
        
        if (curoutput["z"][-1,0,0]<10000) and curoutput["p"][-1,0,0]<10000: # can't have 100hPa at 10km height
            print("WARNING!  Assuming that cesm z has been erroneously divided by 9.8")
            curoutput["z"]*=9.8
        
        outputdata.append(curoutput)
    
    return outputdata

def load_atm(time, info):
    """Load atmospheric variable from a GRIB file"""
    uvfile, scfile = find_atm_file(time, info)
    uvnc_file = grib2nc(uvfile, atmuvlist, info.nc_file_dir)
    scnc_file = grib2nc(scfile, atmvarlist, info.nc_file_dir)

    outputdata = Bunch()
    for s, v in zip(icar_uv_var, atmuvlist):
        nc_data = mygis.read_nc(uvnc_file, v, returnNCvar=True)
        if len(nc_data.data.shape) == 3:
            outputdata[s] = nc_data.data[:, info.ymin : info.ymax, info.xmin : info.xmax]
        else:
            outputdata[s] = nc_data.data[info.ymin : info.ymax, info.xmin : info.xmax]
        nc_data.ncfile.close()

    for s, v in zip(icar_atm_var, atmvarlist):
        nc_data = mygis.read_nc(scnc_file, v, returnNCvar=True)
        if len(nc_data.data.shape) == 3:
            outputdata[s] = nc_data.data[:, info.ymin : info.ymax, info.xmin : info.xmax]
        elif len(nc_data.data.shape) == 2:
            outputdata[s] = nc_data.data[info.ymin : info.ymax, info.xmin : info.xmax]
        elif len(nc_data.data.shape) == 1:
            outputdata[s] = nc_data.data[:]
        else:
            outputdata[s] = nc_data.data.get_value()

        nc_data.ncfile.close()

    return outputdata

def main():
    """docstring for main"""
    plt.figure(figsize=(25,25))
    months=np.arange(1,13)
    for i in range(1,30):
        print(i)
        currentdata=(mygis.read_nc(current_precc_file.format(i+1),"PRECC",returnNCvar=True).data[current_start_point:current_end_point,130:145,200:208].mean(axis=1).mean(axis=1)
                    +mygis.read_nc(current_precl_file.format(i+1),"PRECL",returnNCvar=True).data[current_start_point:current_end_point,130:145,200:208].mean(axis=1).mean(axis=1))
        futuredata =(mygis.read_nc(future_precc_file.format(i+1), "PRECC",returnNCvar=True).data[future_start_point:future_end_point,  130:145,200:208].mean(axis=1).mean(axis=1)
                    +mygis.read_nc(future_precl_file.format(i+1), "PRECL",returnNCvar=True).data[future_start_point:future_end_point,  130:145,200:208].mean(axis=1).mean(axis=1))
        
        cur=np.reshape(currentdata,(10,12)).mean(axis=0)
        fut=np.reshape(futuredata,(10,12)).mean(axis=0)
        delta=1000.*86400.*(fut-cur)
        
        plt.subplot(6,5,i+1)
        # plt.plot(months,delta,label="{:03}".format(i+1),color="black")
        for j in range(12):
            if delta[j]>0:
                plt.fill_between([months[j]-0.4,months[j]+0.4],[delta[j],delta[j]],color="blue")
            else:
                plt.fill_between([months[j]-0.4,months[j]+0.4],[delta[j],delta[j]],color="red")
        plt.plot([0,13],[0,0],":",color="black")
        plt.ylim(-0.5,0.5)
        plt.xlim(0.1,12.9)
        plt.title("Ens {:03}".format(i+1))
        
    # plt.legend()
    plt.tight_layout
    plt.savefig("all_diffs.png")

def vcoord(filename):
    """compute the vertical coordinate in space and time for a given file"""
    na=np.newaxis
    ap= mygis.read_nc(filename,"ap").data[na,:,na,na]
    b = mygis.read_nc(filename,"b").data[na,:,na,na]
    ps= mygis.read_nc(filename,"ps").data[:,na,:,:]
    p= ap+b*ps
    return p

def vcoord(filename):
    """compute the vertical coordinate in space and time for a given file"""
    na=np.newaxis
    a = mygis.read_nc(filename,"lev").data[na,:,na,na]
    b = mygis.read_nc(filename,"b").data[na,:,na,na]
    orog= mygis.read_nc(filename,"orog").data
    z= a+b*orog
    return z

def vcoord(filename):
    """compute the vertical coordinate in space and time for a given file"""
    na=np.newaxis
    ptop = mygis.read_nc(filename,"ptop").data
    sigma = mygis.read_nc(filename,"lev").data[na,:,na,na]
    ps= mygis.read_nc(filename,"ps").data[:,na,:,:]
    p= ptop+sigma*(ps-ptop)
    return p

def main():
    """convert probability coefficients to precip amounts"""
    print("Reading Time data")
    timeseconds=mygis.read_nc(probfile,"time").data
    time=[base_date+datetime.timedelta(i/86400.0) for i in timeseconds]
    print("Reading Spatial data")
    lat=mygis.read_nc(probfile,"latitude").data
    lon=mygis.read_nc(probfile,"longitude").data
    print("Reading Probability regressions")
    prob=mygis.read_nc(probfile,"coefficient").data
    print("Reading Precip regressions")
    prcp=mygis.read_nc(prcpfile,"coefficient").data
    nx=prob.shape[-1]
    ny=prob.shape[-2]
    nt=prob.shape[-3]
    nv=prob.shape[0]
    nt=len(time)
    # print(nv)
    # print(file_variable,variable_name)
    # nt=30
    print(time[0],time[-1],len(time),nt)
    output_data=np.zeros((nt,ny,nx))
    output_data2=np.zeros((nt,ny,nx))
    output_data3=np.zeros((nt,ny,nx))
    output_data4=np.zeros((nt,ny,nx))
    print("starting")
    nvars=(nv-3)/2 # minus 3 for constant, output, residual columns
                   # /2 because we get both the coefficient and the value used for each
    
    for i in range(nt):
        print("  ",i," / ",nt,end="\r")
        sys.stdout.flush()
        # gefs_data=load_gefs(time[i],geo_file=probfile)
        rand_data=load_rand(i)
        curprec=prcp[nvars+1,i,...]
        curprob=prob[nvars+1,i,...]
        # for v in range(1,nv):
        #     # gefs_data[v-1]=gefs_data[v-1]*gains[v-1]+offsets[v-1]
        #     curprec+=gefs_data[v-1]*prcp[v,i,...]
        #     curprob+=gefs_data[v-1]*prob[v,i,...]
        
        # curprec[norm.cdf(rand_data)<curprob]=0
        # output_data[i,...]=curprec
        # output_data2[i,...]=curprob
        curprob/=200.0
        curprob[curprob<norm.cdf(rand_data)]=0
        output_data3[i,...]=curprob
        # anywhere the probability of precip is less than the random number, set precip to 0
        curprec[curprob<norm.cdf(rand_data)]=0
        # add a random component to the precipitation rescaled by the residuals
        curprec+=prcp[nvars+2,i,...] * rand_data
        curprec[curprec<0]=0
        output_data4[i,...]=curprec**(1/3.0)
    
    # mygis.write(output_file,output_data)
    # mygis.write(output_file+"prob",output_data2)
    mygis.write(output_file+"prob_thresh",output_data3)
    mygis.write(output_file+"prec_thresh",output_data4)