Changeset 1992


Ignore:
Timestamp:
Nov 2, 2005, 10:20:13 AM (18 years ago)
Author:
duncan
Message:

adding code to convert gippsland to sww

Location:
inundation/pyvolution
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/data_manager.py

    r1985 r1992  
    4949"""
    5050
    51 from Numeric import concatenate
     51from Numeric import concatenate, array, Float
    5252
    5353from coordinate_transforms.geo_reference import Geo_reference
     
    28392839    outfile.close()
    28402840
     2841
     2842   
     2843def sww2asc_obsolete(basename_in, basename_out = None,
     2844            quantity = None,
     2845            timestep = None,
     2846            reduction = None,
     2847            cellsize = 10,
     2848            verbose = False,
     2849            origin = None):
     2850    """Read SWW file and convert to Digitial Elevation model format (.asc)
     2851
     2852    Example:
     2853
     2854    ncols         3121
     2855    nrows         1800
     2856    xllcorner     722000
     2857    yllcorner     5893000
     2858    cellsize      25
     2859    NODATA_value  -9999
     2860    138.3698 137.4194 136.5062 135.5558 ..........
     2861
     2862    Also write accompanying file with same basename_in but extension .prj
     2863    used to fix the UTM zone, datum, false northings and eastings.
     2864
     2865    The prj format is assumed to be as
     2866
     2867    Projection    UTM
     2868    Zone          56
     2869    Datum         WGS84
     2870    Zunits        NO
     2871    Units         METERS
     2872    Spheroid      WGS84
     2873    Xshift        0.0000000000
     2874    Yshift        10000000.0000000000
     2875    Parameters
     2876
     2877
     2878    if quantity is given, out values from quantity otherwise default to
     2879    elevation
     2880
     2881    if timestep (an index) is given, output quantity at that timestep
     2882
     2883    if reduction is given use that to reduce quantity over all timesteps.
     2884
     2885    """
     2886    from Numeric import array, Float, concatenate, NewAxis, zeros,\
     2887         sometrue
     2888    from utilities.polygon import inside_polygon
     2889
     2890    #FIXME: Should be variable
     2891    datum = 'WGS84'
     2892    false_easting = 500000
     2893    false_northing = 10000000
     2894
     2895    if quantity is None:
     2896        quantity = 'elevation'
     2897
     2898    if reduction is None:
     2899        reduction = max
     2900
     2901    if basename_out is None:
     2902        basename_out = basename_in + '_%s' %quantity
     2903
     2904    swwfile = basename_in + '.sww'
     2905    ascfile = basename_out + '.asc'
     2906    prjfile = basename_out + '.prj'
     2907
     2908
     2909    if verbose: print 'Reading from %s' %swwfile
     2910    #Read sww file
     2911    from Scientific.IO.NetCDF import NetCDFFile
     2912    fid = NetCDFFile(swwfile)
     2913
     2914    #Get extent and reference
     2915    x = fid.variables['x'][:]
     2916    y = fid.variables['y'][:]
     2917    volumes = fid.variables['volumes'][:]
     2918
     2919    ymin = min(y); ymax = max(y)
     2920    xmin = min(x); xmax = max(x)
     2921
     2922    number_of_timesteps = fid.dimensions['number_of_timesteps']
     2923    number_of_points = fid.dimensions['number_of_points']
     2924    if origin is None:
     2925
     2926        #Get geo_reference
     2927        #sww files don't have to have a geo_ref
     2928        try:
     2929            geo_reference = Geo_reference(NetCDFObject=fid)
     2930        except AttributeError, e:
     2931            geo_reference = Geo_reference() #Default georef object
     2932
     2933        xllcorner = geo_reference.get_xllcorner()
     2934        yllcorner = geo_reference.get_yllcorner()
     2935        zone = geo_reference.get_zone()
     2936    else:
     2937        zone = origin[0]
     2938        xllcorner = origin[1]
     2939        yllcorner = origin[2]
     2940
     2941
     2942    #Get quantity and reduce if applicable
     2943    if verbose: print 'Reading quantity %s' %quantity
     2944
     2945    if quantity.lower() == 'depth':
     2946        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
     2947    else:
     2948        q = fid.variables[quantity][:]
     2949
     2950
     2951    if len(q.shape) == 2:
     2952        if verbose: print 'Reducing quantity %s' %quantity
     2953        q_reduced = zeros( number_of_points, Float )
     2954
     2955        for k in range(number_of_points):
     2956            q_reduced[k] = reduction( q[:,k] )
     2957
     2958        q = q_reduced
     2959
     2960    #Now q has dimension: number_of_points
     2961
     2962    #Create grid and update xll/yll corner and x,y
     2963    if verbose: print 'Creating grid'
     2964    ncols = int((xmax-xmin)/cellsize)+1
     2965    nrows = int((ymax-ymin)/cellsize)+1
     2966
     2967    newxllcorner = xmin+xllcorner
     2968    newyllcorner = ymin+yllcorner
     2969
     2970    x = x+xllcorner-newxllcorner
     2971    y = y+yllcorner-newyllcorner
     2972
     2973    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
     2974    assert len(vertex_points.shape) == 2
     2975
     2976
     2977    from Numeric import zeros, Float
     2978    grid_points = zeros ( (ncols*nrows, 2), Float )
     2979
     2980
     2981    for i in xrange(nrows):
     2982        yg = i*cellsize
     2983        for j in xrange(ncols):
     2984            xg = j*cellsize
     2985            k = i*ncols + j
     2986
     2987            grid_points[k,0] = xg
     2988            grid_points[k,1] = yg
     2989
     2990    #Interpolate
     2991    from least_squares import Interpolation
     2992   
     2993
     2994    #FIXME: This should be done with precrop = True, otherwise it'll
     2995    #take forever. With expand_search  set to False, some grid points might
     2996    #miss out....
     2997
     2998    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
     2999                           precrop = False, expand_search = True,
     3000                           verbose = verbose)
     3001
     3002    #Interpolate using quantity values
     3003    if verbose: print 'Interpolating'
     3004    grid_values = interp.interpolate(q).flat
     3005
     3006    #Write
     3007    #Write prj file
     3008    if verbose: print 'Writing %s' %prjfile
     3009    prjid = open(prjfile, 'w')
     3010    prjid.write('Projection    %s\n' %'UTM')
     3011    prjid.write('Zone          %d\n' %zone)
     3012    prjid.write('Datum         %s\n' %datum)
     3013    prjid.write('Zunits        NO\n')
     3014    prjid.write('Units         METERS\n')
     3015    prjid.write('Spheroid      %s\n' %datum)
     3016    prjid.write('Xshift        %d\n' %false_easting)
     3017    prjid.write('Yshift        %d\n' %false_northing)
     3018    prjid.write('Parameters\n')
     3019    prjid.close()
     3020
     3021
     3022   
     3023    if verbose: print 'Writing %s' %ascfile
     3024    NODATA_value = -9999
     3025 
     3026    ascid = open(ascfile, 'w')
     3027
     3028    ascid.write('ncols         %d\n' %ncols)
     3029    ascid.write('nrows         %d\n' %nrows)
     3030    ascid.write('xllcorner     %d\n' %newxllcorner)
     3031    ascid.write('yllcorner     %d\n' %newyllcorner)
     3032    ascid.write('cellsize      %f\n' %cellsize)
     3033    ascid.write('NODATA_value  %d\n' %NODATA_value)
     3034
     3035
     3036    #Get bounding polygon from mesh
     3037    P = interp.mesh.get_boundary_polygon()
     3038    inside_indices = inside_polygon(grid_points, P)
     3039
     3040    for i in range(nrows):
     3041        if verbose and i%((nrows+10)/10)==0:
     3042            print 'Doing row %d of %d' %(i, nrows)
     3043
     3044        for j in range(ncols):
     3045            index = (nrows-i-1)*ncols+j
     3046
     3047            if sometrue(inside_indices == index):
     3048                ascid.write('%f ' %grid_values[index])
     3049            else:
     3050                ascid.write('%d ' %NODATA_value)
     3051
     3052        ascid.write('\n')
     3053
     3054    #Close
     3055    ascid.close()
     3056    fid.close()
     3057
     3058def asc_csiro2sww(bath_file_in, verbose=False):
     3059    pass
     3060
     3061def _read_asc(filename, verbose=False):
     3062    """Read esri file from the following ASCII format (.asc)
     3063
     3064    Example:
     3065
     3066    ncols         3121
     3067    nrows         1800
     3068    xllcorner     722000
     3069    yllcorner     5893000
     3070    cellsize      25
     3071    NODATA_value  -9999
     3072    138.3698 137.4194 136.5062 135.5558 ..........
     3073    """
     3074
     3075    datafile = open(filename)
     3076
     3077    if verbose: print 'Reading DEM from %s' %(filename)
     3078    lines = datafile.readlines()
     3079    datafile.close()
     3080
     3081    if verbose: print 'Got', len(lines), ' lines'
     3082
     3083    ncols = int(lines.pop(0).split()[1].strip())
     3084    nrows = int(lines.pop(0).split()[1].strip())
     3085    xllcorner = float(lines.pop(0).split()[1].strip())
     3086    yllcorner = float(lines.pop(0).split()[1].strip())
     3087    cellsize = float(lines.pop(0).split()[1].strip())
     3088    NODATA_value = int(lines.pop(0).split()[1].strip())
     3089   
     3090    assert len(lines) == nrows
     3091   
     3092    #Store data
     3093    grid = []
     3094   
     3095    n = len(lines)
     3096    for i, line in enumerate(lines):
     3097        cells = line.split()
     3098        assert len(cells) == ncols
     3099        grid.append(array([float(x) for x in cells]))
     3100    grid = array(grid)
     3101
     3102    return xllcorner, yllcorner, cellsize, NODATA_value, grid
     3103   
  • inundation/pyvolution/test_data_manager.py

    r1919 r1992  
    66from Numeric import zeros, array, allclose, Float
    77from util import mean
     8import tempfile
    89
    910from data_manager import *
     
    26532654                verbose = True)
    26542655
    2655 
    2656 
     2656    def test_read_asc(self):
     2657        """Test conversion from dem in ascii format to native NetCDF xya format
     2658        """
     2659
     2660        import time, os
     2661        from Numeric import array, zeros, allclose, Float, concatenate
     2662        from Scientific.IO.NetCDF import NetCDFFile
     2663
     2664        import data_manager
     2665        #Write test asc file
     2666        filename = tempfile.mktemp(".000")
     2667        fid = open(filename, 'w')
     2668        fid.write("""ncols         7
     2669nrows         4
     2670xllcorner     2000.5
     2671yllcorner     3000.5
     2672cellsize      25
     2673NODATA_value  -9999
     2674    97.921    99.285   125.588   180.830   258.645   342.872   415.836
     2675   473.157   514.391   553.893   607.120   678.125   777.283   883.038
     2676   984.494  1040.349  1008.161   900.738   730.882   581.430   514.980
     2677   502.645   516.230   504.739   450.604   388.500   338.097   514.980
     2678""")
     2679        fid.close()
     2680        xllcorner, yllcorner, cellsize, NODATA_value, grid = \
     2681                   data_manager._read_asc(filename, verbose=False)
     2682        self.failUnless(xllcorner  == 2000.5,  'Failed')
     2683        self.failUnless(yllcorner  == 3000.5,  'Failed')
     2684        self.failUnless(cellsize  == 25,  'Failed')
     2685        self.failUnless(NODATA_value  == -9999,  'Failed')
     2686        self.failUnless(grid[0][0]  == 97.921,  'Failed')
     2687        self.failUnless(grid[3][6]  == 514.980,  'Failed')
     2688       
    26572689
    26582690#-------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.