Changeset 5462


Ignore:
Timestamp:
Jul 4, 2008, 11:48:37 AM (14 years ago)
Author:
ole
Message:

John and Ole implemented urs2sts using an ordering file creating smaller and pre-ordered sts files for use with boundary conditions.

Location:
anuga_core/source/anuga/shallow_water
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r5452 r5462  
    49034903def urs2sts(basename_in, basename_out = None, weights=None,
    49044904            verbose = False, origin = None,
    4905             mean_stage=0.0,zscale=1.0,
     4905            mean_stage=0.0, zscale=1.0,
     4906            ordering_filename = None,
    49064907            minlat = None, maxlat = None,
    49074908            minlon = None, maxlon = None):
     
    49134914    origin is a 3-tuple with geo referenced
    49144915    UTM coordinates (zone, easting, northing)
     4916   
     4917    input:
     4918    basename_in: list of source file prefixes
     4919   
     4920    These are combined with the extensions:
     4921    WAVEHEIGHT_MUX2_LABEL = '_waveheight-z-mux2' for stage
     4922    EAST_VELOCITY_MUX2_LABEL =  '_velocity-e-mux2' xmomentum
     4923    NORTH_VELOCITY_MUX2_LABEL =  '_velocity-n-mux2' and ymomentum
     4924   
     4925    to create a 2D list of mux2 file. The rows are associated with each quantity and must have the above extensions
     4926    the columns are the list of file prefixes.
     4927   
     4928    ordering: a .txt file name specifying which mux2 gauge points are sto be stored. This is indicated by the index of
     4929    the gauge in the ordering file.
     4930   
     4931    ordering file format:
     4932    header
     4933    index,longitude,latitude
     4934   
     4935    header='index,longitude,latitude\n'
     4936    If ordering is None all points are taken in the order they appear in the mux2 file subject to rectangular clipping box (see below).
     4937   
     4938   
     4939    minlat, minlon, maxlat, maxlon define a rectangular clipping box and can be used to extract gauge information in a rectangular box. This is useful for extracting
     4940    information at gauge points not on the boundary.
     4941    if ordering_filename is specified clipping box cannot be used
     4942   
     4943   
     4944    output:
     4945    baename_out: name of sts file in which mux2 data is stored.
     4946   
     4947   
    49154948    """
    49164949    import os
     
    49444977    precision = Float
    49454978
     4979    if ordering_filename is not None:
     4980        msg = 'If ordering_filename is specified,'
     4981        msg += ' rectangular clipping box cannot be specified as well. '
     4982        assert minlat is None and maxlat is None and\
     4983            minlon is None and maxlon is None, msg
     4984   
     4985   
    49464986    msg = 'Must use latitudes and longitudes for minlat, maxlon etc'
    4947 
    49484987    if minlat != None:
    49494988        assert -90 < minlat < 90 , msg
     
    49594998            assert maxlon > minlon
    49604999
     5000
     5001    # Check output filename   
    49615002    if basename_out is None:
    4962         stsname = basename_in[0] + '.sts'
    4963     else:
    4964         stsname = basename_out + '.sts'
     5003        msg = 'STS filename must be specified'
     5004        raise Exception, msg
     5005   
     5006    stsname = basename_out + '.sts'
    49655007
    49665008    files_in=[[],[],[]]
     
    49915033    mux={}
    49925034    for i,quantity in enumerate(quantities):
    4993         times_urs, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights)
     5035        # For each quantity read the associated list of source mux2 file with extenstion associated with that quantity
     5036        times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights)
    49945037        if quantity!=quantities[0]:
    49955038            msg='%s, %s and %s have inconsitent gauge data'%(files_in[0],files_in[1],files_in[2])
    4996             assert allclose(times_urs,times_urs_old),msg
     5039            assert allclose(times,times_old),msg
    49975040            assert allclose(latitudes_urs,latitudes_urs_old),msg
    49985041            assert allclose(longitudes_urs,longitudes_urs_old),msg
    49995042            assert allclose(elevation,elevation_old),msg
    50005043            assert allclose(starttime,starttime_old)
    5001         times_urs_old=times_urs
     5044        times_old=times
    50025045        latitudes_urs_old=latitudes_urs
    50035046        longitudes_urs_old=longitudes_urs
    50045047        elevation_old=elevation
    50055048        starttime_old=starttime
    5006        
    5007     if (minlat is not None) and (minlon is not None) and (maxlat is not None) and (maxlon is not None):
     5049
     5050       
     5051    if ordering_filename is not None:
     5052        # Read ordering file
     5053       
     5054        try:
     5055            fid=open(ordering_filename, 'r')
     5056            file_header=fid.readline().split(',')
     5057            ordering_lines=fid.readlines()
     5058            fid.close()
     5059        except:
     5060            msg = 'Cannot open %s'%ordering_filename
     5061            raise Exception, msg
     5062
     5063        reference_header = 'index, longitude, latitude\n'.split(',')
     5064        for i in range(3):
     5065            if not file_header[i] == reference_header[i]:
     5066                msg = 'File must contain header\n'+header+'\n'
     5067                raise Exception, msg
     5068
     5069        if len(ordering_lines)<2:
     5070            msg = 'File must contain at least two points'
     5071            raise Exception, msg
     5072
     5073   
     5074       
     5075        permutation = [int(line.split(',')[0]) for line in ordering_lines]
     5076   
     5077        latitudes = take(latitudes_urs, permutation)
     5078        longitudes = take(longitudes_urs, permutation)
     5079       
     5080        # Self check - can be removed to improve speed
     5081        ref_longitudes = [float(line.split(',')[1]) for line in ordering_lines]               
     5082        ref_latitudes = [float(line.split(',')[2]) for line in ordering_lines]       
     5083       
     5084        msg = 'Longitudes specified in ordering file do not match those found in mux files'
     5085        msg += 'I got %s instead of %s (only beginning shown)'\
     5086            %(str(longitudes[:10]) + '...',
     5087              str(ref_longitudes[:10]) + '...')       
     5088        assert allclose(longitudes, ref_longitudes), msg       
     5089       
     5090        msg = 'Latitudes specified in ordering file do not match those found in mux files. '
     5091        msg += 'I got %s instead of %s (only beginning shown)'\
     5092            %(str(latitudes[:10]) + '...',
     5093              str(ref_latitudes[:10]) + '...')               
     5094        assert allclose(latitudes, ref_latitudes), msg
     5095       
     5096    elif (minlat is not None) and (minlon is not None) and (maxlat is not None) and (maxlon is not None):
     5097        #FIXME - this branch is probably not working   
    50085098        if verbose: print 'Cliping urs data'
    50095099        latitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),latitudes_urs)
    50105100        longitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),longitudes_urs)
    5011         times = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),times_urs)
     5101       
     5102        permutation = range(latitudes.shape[0])       
     5103       
     5104        msg = 'Clipping is not done yet'
     5105        raise Exception, msg
     5106       
    50125107    else:
    50135108        latitudes=latitudes_urs
    50145109        longitudes=longitudes_urs
    5015         times=times_urs
     5110        permutation = range(latitudes.shape[0])               
     5111       
    50165112
    50175113    msg='File is empty and or clipped region not in file region'
     
    50325128    #starttime = times[0]
    50335129    sts = Write_sts()
    5034     sts.store_header(outfile, times+starttime,number_of_points, description=description,
    5035                      verbose=verbose,sts_precision=Float)
     5130    sts.store_header(outfile,
     5131                     times+starttime,
     5132                     number_of_points,
     5133                     description=description,
     5134                     verbose=verbose,
     5135                     sts_precision=Float)
    50365136   
    50375137    # Store
     
    50675167    stage = outfile.variables['stage']
    50685168    xmomentum = outfile.variables['xmomentum']
    5069     ymomentum =outfile.variables['ymomentum']
     5169    ymomentum = outfile.variables['ymomentum']
    50705170
    50715171    if verbose: print 'Converting quantities'
    50725172    for j in range(len(times)):
    5073         for i in range(number_of_points):
    5074             w = zscale*mux['HA'][i,j] + mean_stage
    5075             h=w-elevation[i]
     5173        for i, index in enumerate(permutation):
     5174            w = zscale*mux['HA'][index,j] + mean_stage
     5175            h=w-elevation[index]
    50765176            stage[j,i] = w
    5077             #delete following two lines once velcotu files are read in.
    5078             xmomentum[j,i] = mux['UA'][i,j]*h
    5079             ymomentum[j,i] = mux['VA'][i,j]*h
     5177
     5178            xmomentum[j,i] = mux['UA'][index,j]*h
     5179            ymomentum[j,i] = mux['VA'][index,j]*h
    50805180
    50815181    outfile.close()
     
    51195219    xllcorner = fid.xllcorner[0]
    51205220    yllcorner = fid.yllcorner[0]
    5121     #Points stored in sts file are normailsed to [xllcorner,yllcorner] but
     5221    #Points stored in sts file are normalised to [xllcorner,yllcorner] but
    51225222    #we cannot assume that boundary polygon will be. At least the
    51235223    #additional points specified by the user after this function is called
     
    51325232    yllcorner = fid.yllcorner[0]
    51335233
     5234    # Walk through ordering file
    51345235    boundary_polygon=[]
    51355236    for line in lines:
     
    51395240            zone,easting,northing=redfearn(float(fields[1]),float(fields[2]))
    51405241            if not zone == fid.zone[0]:
    5141                 msg = 'Inconsitent zones'
     5242                msg = 'Inconsistent zones'
    51425243                raise Exception, msg
    51435244        else:
    51445245            easting = float(fields[1])
    51455246            northing = float(fields[2])
     5247           
    51465248        if not allclose(array([easting,northing]),sts_points[index]):
    51475249            msg = "Points in sts file do not match those in the"+\
    5148                 ".txt file spcifying the order of the boundary points"
     5250                ".txt file specifying the order of the boundary points"
    51495251            raise Exception, msg
    51505252        boundary_polygon.append(sts_points[index].tolist())
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5442 r5462  
    63196319                                      va=va)
    63206320
    6321         urs2sts(base_name,mean_stage=tide,verbose=False)
     6321        urs2sts(base_name,
     6322                basename_out=base_name,
     6323                mean_stage=tide,verbose=False)
    63226324
    63236325        # now I want to check the sts file ...
     
    64416443
    64426444        # Call urs2sts with multiple mux files
    6443         urs2sts([base_nameI, base_nameII], weights=[1.0,1.0],
     6445        urs2sts([base_nameI, base_nameII],
     6446                basename_out=base_nameI,
     6447                weights=[1.0,1.0],
    64446448                mean_stage=tide,
    64456449                verbose=False)
     
    65346538        self.delete_mux(filesII)
    65356539        os.remove(sts_file)
     6540       
     6541       
     6542    def test_urs2sts_ordering(self):
     6543        """Test multiple sources with ordering file
     6544        """
     6545       
     6546        tide=0
     6547        time_step_count = 3
     6548        time_step = 2
     6549        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
     6550        n=len(lat_long_points)
     6551        first_tstep=ones(n,Int)
     6552        first_tstep[0]+=1
     6553        first_tstep[2]+=1
     6554        last_tstep=(time_step_count)*ones(n,Int)
     6555        last_tstep[0]-=1
     6556
     6557        gauge_depth=20*ones(n,Float)
     6558        ha=2*ones((n,time_step_count),Float)
     6559        ha[0]=arange(0,time_step_count)
     6560        ha[1]=arange(time_step_count,2*time_step_count)
     6561        ha[2]=arange(2*time_step_count,3*time_step_count)
     6562        ha[3]=arange(3*time_step_count,4*time_step_count)
     6563        ua=5*ones((n,time_step_count),Float)
     6564        va=-10*ones((n,time_step_count),Float)
     6565
     6566        # Create two identical mux files to be combined by urs2sts
     6567        base_nameI, filesI = self.write_mux2(lat_long_points,
     6568                                             time_step_count, time_step,
     6569                                             first_tstep, last_tstep,
     6570                                             depth=gauge_depth,
     6571                                             ha=ha,
     6572                                             ua=ua,
     6573                                             va=va)
     6574
     6575        base_nameII, filesII = self.write_mux2(lat_long_points,
     6576                                               time_step_count, time_step,
     6577                                               first_tstep, last_tstep,
     6578                                               depth=gauge_depth,
     6579                                               ha=ha,
     6580                                               ua=ua,
     6581                                               va=va)
     6582
     6583                                               
     6584        # Create ordering file
     6585        permutation = [3,0,2]
     6586
     6587        # Do it wrongly at first and check that exception is being raised
     6588        _, ordering_filename = tempfile.mkstemp('')
     6589        order_fid = open(ordering_filename, 'w') 
     6590        order_fid.write('index, longitude, latitude\n')
     6591        for index in permutation:
     6592            order_fid.write('%d, %f, %f\n' %(index,
     6593                                             lat_long_points[index][0],
     6594                                             lat_long_points[index][1]))
     6595        order_fid.close()
     6596       
     6597        try:
     6598            urs2sts([base_nameI, base_nameII],
     6599                    basename_out=base_nameI,
     6600                    ordering_filename=ordering_filename,
     6601                    weights=[1.0,1.0],
     6602                    mean_stage=tide,
     6603                    verbose=False) 
     6604            os.remove(ordering_filename)           
     6605        except:
     6606            pass
     6607        else:
     6608            msg = 'Should have caught wrong lat longs'
     6609            raise Exception, msg
     6610
     6611
     6612       
     6613           
     6614        # Then do it right
     6615        _, ordering_filename = tempfile.mkstemp('')
     6616        order_fid = open(ordering_filename, 'w') 
     6617        order_fid.write('index, longitude, latitude\n')
     6618        for index in permutation:
     6619            order_fid.write('%d, %f, %f\n' %(index,
     6620                                             lat_long_points[index][1],
     6621                                             lat_long_points[index][0]))
     6622        order_fid.close()
     6623       
     6624           
     6625
     6626                                               
     6627        # Call urs2sts with multiple mux files
     6628        urs2sts([base_nameI, base_nameII],
     6629                basename_out=base_nameI,
     6630                ordering_filename=ordering_filename,
     6631                weights=[1.0,1.0],
     6632                mean_stage=tide,
     6633                verbose=False)
     6634
     6635        # now I want to check the sts file ...
     6636        sts_file = base_nameI + '.sts'
     6637
     6638        #Let's interrogate the sts file
     6639        # Note, the sts info is not gridded.  It is point data.
     6640        fid = NetCDFFile(sts_file)
     6641
     6642        # Make x and y absolute
     6643        x = fid.variables['x'][:]
     6644        y = fid.variables['y'][:]
     6645
     6646        geo_reference = Geo_reference(NetCDFObject=fid)
     6647        points = geo_reference.get_absolute(map(None, x, y))
     6648        points = ensure_numeric(points)
     6649
     6650        x = points[:,0]
     6651        y = points[:,1]
     6652
     6653        for i, index in enumerate(permutation):
     6654            # Check that STS points are stored in the correct order
     6655           
     6656            # Work out the UTM coordinates sts point i
     6657            zone, e, n = redfearn(lat_long_points[index][0],
     6658                                  lat_long_points[index][1])             
     6659           
     6660            assert allclose([x[i],y[i]], [e,n])
     6661           
     6662                       
     6663        # Check the time vector
     6664        times = fid.variables['time'][:]
     6665
     6666        times_actual = []
     6667        for i in range(time_step_count):
     6668            times_actual.append(time_step * i)
     6669
     6670        assert allclose(ensure_numeric(times),
     6671                        ensure_numeric(times_actual))
     6672                       
     6673
     6674        # Check sts values
     6675        stage = fid.variables['stage'][:]
     6676        xmomentum = fid.variables['xmomentum'][:]
     6677        ymomentum = fid.variables['ymomentum'][:]
     6678        elevation = fid.variables['elevation'][:]
     6679
     6680        # Set original data used to write mux file to be zero when gauges are
     6681        # not recdoring
     6682       
     6683        ha[0][0]=0.0
     6684        ha[0][time_step_count-1]=0.0
     6685        ha[2][0]=0.0
     6686        ua[0][0]=0.0
     6687        ua[0][time_step_count-1]=0.0
     6688        ua[2][0]=0.0
     6689        va[0][0]=0.0
     6690        va[0][time_step_count-1]=0.0
     6691        va[2][0]=0.0;
     6692
     6693       
     6694        # The stage stored in the .sts file should be the sum of the stage
     6695        # in the two mux2 files because both have weights = 1. In this case
     6696        # the mux2 files are the same so stage == 2.0 * ha
     6697        #print 2.0*transpose(ha) - stage
     6698       
     6699        ha_permutation = take(ha, permutation)
     6700        ua_permutation = take(ua, permutation)         
     6701        va_permutation = take(va, permutation)                 
     6702        gauge_depth_permutation = take(gauge_depth, permutation)                         
     6703
     6704       
     6705        assert allclose(2.0*transpose(ha_permutation), stage)  # Meters
     6706
     6707        #Check the momentums - ua
     6708        #momentum = velocity*(stage-elevation)
     6709        # elevation = - depth
     6710        #momentum = velocity_ua *(stage+depth)
     6711
     6712        depth=zeros((len(lat_long_points),time_step_count),Float)
     6713        for i in range(len(lat_long_points)):
     6714            depth[i]=gauge_depth[i]+tide+2.0*ha[i]
     6715            #2.0*ha necessary because using two files with weights=1 are used
     6716           
     6717        depth_permutation = take(depth, permutation)                     
     6718       
     6719
     6720        # The xmomentum stored in the .sts file should be the sum of the ua
     6721        # in the two mux2 files multiplied by the depth.
     6722        assert allclose(2.0*transpose(ua_permutation*depth_permutation), xmomentum)
     6723
     6724        #Check the momentums - va
     6725        #momentum = velocity*(stage-elevation)
     6726        # elevation = - depth
     6727        #momentum = velocity_va *(stage+depth)
     6728
     6729        # The ymomentum stored in the .sts file should be the sum of the va
     6730        # in the two mux2 files multiplied by the depth.
     6731        assert allclose(2.0*transpose(va_permutation*depth_permutation), ymomentum)
     6732
     6733        # check the elevation values.
     6734        # -ve since urs measures depth, sww meshers height,
     6735        assert allclose(-gauge_depth_permutation, elevation)  #Meters
     6736
     6737        fid.close()
     6738        self.delete_mux(filesI)
     6739        self.delete_mux(filesII)
     6740        os.remove(sts_file)
     6741       
    65366742
    65376743    def test_file_boundary_stsI(self):
     
    65636769
    65646770        sts_file=base_name
    6565         urs2sts(base_name,sts_file,mean_stage=tide,verbose=False)
     6771        urs2sts(base_name,
     6772                sts_file,
     6773                mean_stage=tide,verbose=False)
    65666774        self.delete_mux(files)
    65676775
Note: See TracChangeset for help on using the changeset viewer.