Changeset 1833


Ignore:
Timestamp:
Sep 15, 2005, 8:53:36 AM (19 years ago)
Author:
tdhu
Message:

Testing is as yet uncomplete. Test_data_manager is currently running data_manager to import an sww to a ermapper grid.

Location:
inundation/pyvolution
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/data_manager.py

    r1797 r1833  
    27092709    outfile.close()
    27102710
     2711def sww2ers(basename_in, basename_out = None,
     2712            quantity = None,
     2713            timestep = None,
     2714            reduction = None,
     2715            cellsize = 10,
     2716            verbose = False,
     2717            origin = None,
     2718            datum = 'WGS84'):
     2719   
     2720    """Read SWW file and convert to Digitial Elevation model format (.asc)
     2721
     2722    Example:
     2723
     2724    ncols         3121
     2725    nrows         1800
     2726    xllcorner     722000
     2727    yllcorner     5893000
     2728    cellsize      25
     2729    NODATA_value  -9999
     2730    138.3698 137.4194 136.5062 135.5558 ..........
     2731
     2732    Also write accompanying file with same basename_in but extension .prj
     2733    used to fix the UTM zone, datum, false northings and eastings.
     2734
     2735    The prj format is assumed to be as
     2736
     2737    Projection    UTM
     2738    Zone          56
     2739    Datum         WGS84
     2740    Zunits        NO
     2741    Units         METERS
     2742    Spheroid      WGS84
     2743    Xshift        0.0000000000
     2744    Yshift        10000000.0000000000
     2745    Parameters
     2746
     2747
     2748    if quantity is given, out values from quantity otherwise default to
     2749    elevation
     2750
     2751    if timestep (an index) is given, output quantity at that timestep
     2752
     2753    if reduction is given use that to reduce quantity over all timesteps.
     2754
     2755    """
     2756    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
     2757    import ermapper_grids
     2758
     2759    header = {}
     2760    false_easting = 500000
     2761    false_northing = 10000000
     2762    NODATA_value = 0
     2763
     2764    if quantity is None:
     2765        quantity = 'elevation'
     2766
     2767    if reduction is None:
     2768        reduction = max
     2769
     2770    if basename_out is None:
     2771        basename_out = basename_in + '_%s' %quantity
     2772 
     2773    swwfile = basename_in + '.sww'
     2774    # Note the use of a .ers extension is optional (write_ermapper_grid will
     2775    # deal with either option
     2776    ersfile = basename_out
     2777##    prjfile = basename_out + '.prj'
     2778
     2779
     2780    if verbose: print 'Reading from %s' %swwfile
     2781    #Read sww file
     2782    from Scientific.IO.NetCDF import NetCDFFile
     2783    fid = NetCDFFile(swwfile)
     2784
     2785    #Get extent and reference
     2786    x = fid.variables['x'][:]
     2787    y = fid.variables['y'][:]
     2788    volumes = fid.variables['volumes'][:]
     2789
     2790    ymin = min(y); ymax = max(y)
     2791    xmin = min(x); xmax = max(x)
     2792
     2793    number_of_timesteps = fid.dimensions['number_of_timesteps']
     2794    number_of_points = fid.dimensions['number_of_points']
     2795    if origin is None:
     2796
     2797        #Get geo_reference
     2798        #sww files don't have to have a geo_ref
     2799        try:
     2800            geo_reference = Geo_reference(NetCDFObject=fid)
     2801        except AttributeError, e:
     2802            geo_reference = Geo_reference() #Default georef object
     2803
     2804        xllcorner = geo_reference.get_xllcorner()
     2805        yllcorner = geo_reference.get_yllcorner()
     2806        zone = geo_reference.get_zone()
     2807    else:
     2808        zone = origin[0]
     2809        xllcorner = origin[1]
     2810        yllcorner = origin[2]
     2811
     2812
     2813    #Get quantity and reduce if applicable
     2814    if verbose: print 'Reading quantity %s' %quantity
     2815
     2816    if quantity.lower() == 'depth':
     2817        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
     2818    else:
     2819        q = fid.variables[quantity][:]
     2820
     2821
     2822    if len(q.shape) == 2:
     2823        if verbose: print 'Reducing quantity %s' %quantity
     2824        q_reduced = zeros( number_of_points, Float )
     2825
     2826        for k in range(number_of_points):
     2827            q_reduced[k] = reduction( q[:,k] )
     2828
     2829        q = q_reduced
     2830
     2831    #Now q has dimension: number_of_points
     2832
     2833    #Create grid and update xll/yll corner and x,y
     2834    if verbose: print 'Creating grid'
     2835    ncols = int((xmax-xmin)/cellsize)+1
     2836    nrows = int((ymax-ymin)/cellsize)+1
     2837
     2838    newxllcorner = xmin+xllcorner
     2839    newyllcorner = ymin+yllcorner
     2840
     2841    x = x+xllcorner-newxllcorner
     2842    y = y+yllcorner-newyllcorner
     2843
     2844    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
     2845    assert len(vertex_points.shape) == 2
     2846
     2847
     2848    from Numeric import zeros, Float
     2849    grid_points = zeros ( (ncols*nrows, 2), Float )
     2850
     2851
     2852    for i in xrange(nrows):
     2853##        yg = i*cellsize
     2854        yg = (nrows-i)*cellsize # this will flip the order of the y values
     2855        for j in xrange(ncols):
     2856            xg = j*cellsize
     2857            k = i*ncols + j
     2858
     2859##            print k, xg, yg
     2860            grid_points[k,0] = xg
     2861            grid_points[k,1] = yg
     2862
     2863    #Interpolate
     2864    from least_squares import Interpolation
     2865    from util import inside_polygon
     2866
     2867    #FIXME: This should be done with precrop = True, otherwise it'll
     2868    #take forever. With expand_search  set to False, some grid points might
     2869    #miss out....
     2870
     2871    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
     2872                           precrop = False, expand_search = True,
     2873                           verbose = verbose)
     2874
     2875    #Interpolate using quantity values
     2876    if verbose: print 'Interpolating'
     2877    grid_values = interp.interpolate(q).flat
     2878    grid_values = reshape(grid_values,(nrows, ncols))
     2879
     2880
     2881    # setup header information
     2882    header['datum'] = '"' + datum + '"'
     2883    # FIXME The use of hardwired UTM and zone number needs to be made optional
     2884    # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
     2885    header['projection'] = '"UTM-' + str(zone) + '"'
     2886    header['coordinatetype'] = 'EN'
     2887    if header['coordinatetype'] == 'LL':
     2888        header['longitude'] = str(newxllcorner)
     2889        header['latitude'] = str(newyllcorner)
     2890    elif header['coordinatetype'] == 'EN':
     2891        header['eastings'] = str(newxllcorner)
     2892        header['northings'] = str(newyllcorner)
     2893    header['nullcellvalue'] = str(NODATA_value)
     2894    header['xdimension'] = str(cellsize)
     2895    header['ydimension'] = str(cellsize)
     2896    header['value'] = '"' + quantity + '"'
     2897
     2898
     2899    #Write
     2900    if verbose: print 'Writing %s' %ersfile
     2901    ermapper_grids.write_ermapper_grid(ersfile,grid_values, header)
     2902
     2903##    ascid = open(ascfile, 'w')
     2904##
     2905##    ascid.write('ncols         %d\n' %ncols)
     2906##    ascid.write('nrows         %d\n' %nrows)
     2907##    ascid.write('xllcorner     %d\n' %newxllcorner)
     2908##    ascid.write('yllcorner     %d\n' %newyllcorner)
     2909##    ascid.write('cellsize      %f\n' %cellsize)
     2910##    ascid.write('NODATA_value  %d\n' %NODATA_value)
     2911##
     2912##
     2913##    #Get bounding polygon from mesh
     2914##    P = interp.mesh.get_boundary_polygon()
     2915##    inside_indices = inside_polygon(grid_points, P)
     2916##
     2917##    for i in range(nrows):
     2918##        if verbose and i%((nrows+10)/10)==0:
     2919##            print 'Doing row %d of %d' %(i, nrows)
     2920##
     2921##        for j in range(ncols):
     2922##            index = (nrows-i-1)*ncols+j
     2923##
     2924##            if sometrue(inside_indices == index):
     2925##                ascid.write('%f ' %grid_values[index])
     2926##            else:
     2927##                ascid.write('%d ' %NODATA_value)
     2928##
     2929##        ascid.write('\n')
     2930##
     2931##    #Close
     2932##    ascid.close()
     2933##    fid.close()
  • inundation/pyvolution/ermapper_grids.py

    r1813 r1833  
    55
    66def write_ermapper_grid(ofile, data, header = {}):
    7     # function write_ermapper_grid(ofile, data, header = {}):
    8     #
    9     # Function to write a 2D Numeric array to an ERMapper grid
    10     #
    11     # Input Parameters:
    12     # ofile
    13 
     7    """
     8    write_ermapper_grid(ofile, data, header = {}):
     9   
     10    Function to write a 2D Numeric array to an ERMapper grid.  There are a series of conventions adopted within
     11    this code, specifically:
     12    1)  The registration coordinate for the data is the SW (or lower-left) corner of the data
     13    2)  The registration coordinates refer to cell centres
     14    3)  The data is a 2D Numeric array with the NW-most data in element (0,0) and the SE-most data in element (N,M)
     15        where N is the last line and M is the last column
     16    4)  There has been no testng of the use of a rotated grid.  Best to keep data in an NS orientation
     17   
     18    Input Parameters:
     19    ofile:      string - filename for output (note the output will consist of two files
     20                ofile and ofile.ers.  Either of these can be entered into this function
     21    data:       Numeric.array - 2D array containing the data to be output to the grid
     22    header:     dictionary - contains spatial information about the grid, in particular:
     23                    header['datum'] datum for the data ('"GDA94"')
     24                    header['projection'] - either '"GEOGRAPHIC"' or '"PROJECTED"'
     25                    header['coordinatetype'] - either 'EN' (for eastings/northings) or
     26                                                      'LL' (for lat/long)
     27                    header['rotation'] - rotation of grid ('0:0:0.0')
     28                    header['celltype'] - data type for writing data ('IEEE4ByteReal')
     29                    header['nullcellvalue'] - value for null cells ('-99999')
     30                    header['xdimension'] - cell size in x-dir in units dictated by 'coordinatetype' ('100')
     31                    header['registrationcellx'] == '0'
     32                    header['ydimension'] - cell size in y-dir in units dictated by 'coordinatetype' ('100')
     33                    header['longitude'] - co-ordinate of registration cell ('0:0:0')
     34                    header['latitude'] - co-ordinate of registration line ('0:0:0')
     35                    header['nrofbands'] - number of bands ('1')
     36                    header['value'] - name of grid ('"Default_Band"')
     37
     38                    Some entries are determined automatically from the data
     39                    header['nroflines'] - number of lines in data
     40                    header['nrofcellsperline'] - number of columns in data
     41                    header['registrationcelly'] == last line of data
     42    """
    1443    # extract filenames for header and data files from ofile
    1544    ers_index = ofile.find('.ers')
     
    79108    fid.write('\t\tCoordinateType\t = ' + header['coordinatetype'] + '\n')
    80109    fid.write('\t\tRotation\t\t = ' + header['rotation'] + '\n')
     110    fid.write('\t\tUnits\t\t = ' + header['units'] + '\n')
    81111    fid.write('\tCoordinateSpace End\n')
    82112
     
    97127    # Write the registration coordinate information
    98128    fid.write('\t\tRegistrationCoord Begin\n')
     129    print X_Class
    99130    fid.write('\t\t\t' + X_Class + '\t\t\t = ' + header[X_Class.lower()] + '\n')
    100131    fid.write('\t\t\t' + Y_Class + '\t\t\t = ' + header[Y_Class.lower()] + '\n')
     
    163194    if not header.has_key('rotation'):
    164195        header['rotation'] = '0:0:0.0'
     196    if not header.has_key('units'):
     197        header['units'] = '"METERS"'
    165198    if not header.has_key('celltype'):
    166199        header['celltype'] = 'IEEE4ByteReal'
  • inundation/pyvolution/test_data_manager.py

    r1797 r1833  
    20522052        os.remove(root + '_100.dem')
    20532053
     2054    def test_sww2ers_simple(self):
     2055        """Test that sww information can be converted correctly to asc/prj
     2056        format readable by e.g. ArcView
     2057        """
     2058
     2059        import time, os
     2060        from Numeric import array, zeros, allclose, Float, concatenate
     2061        from Scientific.IO.NetCDF import NetCDFFile
     2062
     2063        #Setup
     2064        self.domain.filename = 'datatest'
     2065
     2066        headerfile = self.domain.filename + '.ers'
     2067        swwfile = self.domain.filename + '.sww'
     2068
     2069        self.domain.set_datadir('.')
     2070        self.domain.format = 'sww'
     2071        self.domain.smooth = True
     2072        self.domain.set_quantity('elevation', lambda x,y: -x-y)
     2073
     2074        self.domain.geo_reference = Geo_reference(56,308500,6189000)
     2075
     2076        sww = get_dataobject(self.domain)
     2077        sww.store_connectivity()
     2078        sww.store_timestep('stage')
     2079
     2080        self.domain.evolve_to_end(finaltime = 0.01)
     2081        sww.store_timestep('stage')
     2082
     2083        cellsize = 0.25
     2084        #Check contents
     2085        #Get NetCDF
     2086
     2087        fid = NetCDFFile(sww.filename, 'r')
     2088
     2089        # Get the variables
     2090        x = fid.variables['x'][:]
     2091        y = fid.variables['y'][:]
     2092        z = fid.variables['elevation'][:]
     2093        time = fid.variables['time'][:]
     2094        stage = fid.variables['stage'][:]
     2095
     2096
     2097        #Export to ascii/prj files
     2098        sww2ers(self.domain.filename,
     2099                quantity = 'elevation',
     2100                cellsize = cellsize,
     2101                verbose = False)
     2102
     2103
     2104##        #Check prj (meta data)
     2105##        prjid = open(prjfile)
     2106##        lines = prjid.readlines()
     2107##        prjid.close()
     2108##
     2109##        L = lines[0].strip().split()
     2110##        assert L[0].strip().lower() == 'projection'
     2111##        assert L[1].strip().lower() == 'utm'
     2112##
     2113##        L = lines[1].strip().split()
     2114##        assert L[0].strip().lower() == 'zone'
     2115##        assert L[1].strip().lower() == '56'
     2116##
     2117##        L = lines[2].strip().split()
     2118##        assert L[0].strip().lower() == 'datum'
     2119##        assert L[1].strip().lower() == 'wgs84'
     2120##
     2121##        L = lines[3].strip().split()
     2122##        assert L[0].strip().lower() == 'zunits'
     2123##        assert L[1].strip().lower() == 'no'
     2124##
     2125##        L = lines[4].strip().split()
     2126##        assert L[0].strip().lower() == 'units'
     2127##        assert L[1].strip().lower() == 'meters'
     2128##
     2129##        L = lines[5].strip().split()
     2130##        assert L[0].strip().lower() == 'spheroid'
     2131##        assert L[1].strip().lower() == 'wgs84'
     2132##
     2133##        L = lines[6].strip().split()
     2134##        assert L[0].strip().lower() == 'xshift'
     2135##        assert L[1].strip().lower() == '500000'
     2136##
     2137##        L = lines[7].strip().split()
     2138##        assert L[0].strip().lower() == 'yshift'
     2139##        assert L[1].strip().lower() == '10000000'
     2140##
     2141##        L = lines[8].strip().split()
     2142##        assert L[0].strip().lower() == 'parameters'
     2143##
     2144##
     2145##        #Check asc file
     2146##        ascid = open(ascfile)
     2147##        lines = ascid.readlines()
     2148##        ascid.close()
     2149##
     2150##        L = lines[0].strip().split()
     2151##        assert L[0].strip().lower() == 'ncols'
     2152##        assert L[1].strip().lower() == '5'
     2153##
     2154##        L = lines[1].strip().split()
     2155##        assert L[0].strip().lower() == 'nrows'
     2156##        assert L[1].strip().lower() == '5'
     2157##
     2158##        L = lines[2].strip().split()
     2159##        assert L[0].strip().lower() == 'xllcorner'
     2160##        assert allclose(float(L[1].strip().lower()), 308500)
     2161##
     2162##        L = lines[3].strip().split()
     2163##        assert L[0].strip().lower() == 'yllcorner'
     2164##        assert allclose(float(L[1].strip().lower()), 6189000)
     2165##
     2166##        L = lines[4].strip().split()
     2167##        assert L[0].strip().lower() == 'cellsize'
     2168##        assert allclose(float(L[1].strip().lower()), cellsize)
     2169##
     2170##        L = lines[5].strip().split()
     2171##        assert L[0].strip() == 'NODATA_value'
     2172##        assert L[1].strip().lower() == '-9999'
     2173##
     2174##        #Check grid values
     2175##        for j in range(5):
     2176##            L = lines[6+j].strip().split()
     2177##            y = (4-j) * cellsize
     2178##            for i in range(5):
     2179##                assert allclose(float(L[i]), -i*cellsize - y)
     2180##
     2181
     2182        fid.close()
     2183        print fid
     2184        #Cleanup
     2185        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
     2186        #os.remove(sww.filename)
     2187
     2188    def testz_sww2ers_real(self):
     2189        """Test that sww information can be converted correctly to asc/prj
     2190        format readable by e.g. ArcView
     2191        """
     2192
     2193        import time, os
     2194        from Numeric import array, zeros, allclose, Float, concatenate
     2195        from Scientific.IO.NetCDF import NetCDFFile
     2196
     2197
     2198
     2199
     2200        #Export to ascii/prj files
     2201        sww2ers('karratha_100m.sww',
     2202                quantity = 'depth',
     2203                cellsize = 5,
     2204                verbose = False)
     2205
     2206
     2207
    20542208
    20552209#-------------------------------------------------------------
    20562210if __name__ == "__main__":
    2057     suite = unittest.makeSuite(Test_Data_Manager,'test')
     2211    suite = unittest.makeSuite(Test_Data_Manager,'testz')
    20582212    #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
    20592213    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
Note: See TracChangeset for help on using the changeset viewer.