Changeset 1985


Ignore:
Timestamp:
Oct 31, 2005, 10:38:30 PM (18 years ago)
Author:
ole
Message:

Removed code made obsolete 4 October 2005 (changeset:1866)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • inundation/pyvolution/data_manager.py

    r1984 r1985  
    28392839    outfile.close()
    28402840
    2841 
    2842    
    2843 def 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 
    3058 def sww2ers_obsolete(basename_in, basename_out = None,
    3059             quantity = None,
    3060             timestep = None,
    3061             reduction = None,
    3062             cellsize = 10,
    3063             verbose = False,
    3064             origin = None,
    3065             datum = 'WGS84'):
    3066    
    3067     """Read SWW file and convert to Digitial Elevation model format (.asc)
    3068 
    3069     Example:
    3070 
    3071     ncols         3121
    3072     nrows         1800
    3073     xllcorner     722000
    3074     yllcorner     5893000
    3075     cellsize      25
    3076     NODATA_value  -9999
    3077     138.3698 137.4194 136.5062 135.5558 ..........
    3078 
    3079     Also write accompanying file with same basename_in but extension .prj
    3080     used to fix the UTM zone, datum, false northings and eastings.
    3081 
    3082     The prj format is assumed to be as
    3083 
    3084     Projection    UTM
    3085     Zone          56
    3086     Datum         WGS84
    3087     Zunits        NO
    3088     Units         METERS
    3089     Spheroid      WGS84
    3090     Xshift        0.0000000000
    3091     Yshift        10000000.0000000000
    3092     Parameters
    3093 
    3094 
    3095     if quantity is given, out values from quantity otherwise default to
    3096     elevation
    3097 
    3098     if timestep (an index) is given, output quantity at that timestep
    3099 
    3100     if reduction is given use that to reduce quantity over all timesteps.
    3101 
    3102     """
    3103     from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
    3104     import ermapper_grids
    3105     from utilities.polygon import inside_polygon   
    3106 
    3107     header = {}
    3108     false_easting = 500000
    3109     false_northing = 10000000
    3110     NODATA_value = 0
    3111 
    3112     if quantity is None:
    3113         quantity = 'elevation'
    3114 
    3115     if reduction is None:
    3116         reduction = max
    3117 
    3118     if basename_out is None:
    3119         basename_out = basename_in + '_%s' %quantity
    3120  
    3121     swwfile = basename_in + '.sww'
    3122     # Note the use of a .ers extension is optional (write_ermapper_grid will
    3123     # deal with either option
    3124     ersfile = basename_out
    3125 ##    prjfile = basename_out + '.prj'
    3126 
    3127 
    3128     if verbose: print 'Reading from %s' %swwfile
    3129     #Read sww file
    3130     from Scientific.IO.NetCDF import NetCDFFile
    3131     fid = NetCDFFile(swwfile)
    3132 
    3133     #Get extent and reference
    3134     x = fid.variables['x'][:]
    3135     y = fid.variables['y'][:]
    3136     volumes = fid.variables['volumes'][:]
    3137 
    3138     ymin = min(y); ymax = max(y)
    3139     xmin = min(x); xmax = max(x)
    3140 
    3141     number_of_timesteps = fid.dimensions['number_of_timesteps']
    3142     number_of_points = fid.dimensions['number_of_points']
    3143     if origin is None:
    3144 
    3145         #Get geo_reference
    3146         #sww files don't have to have a geo_ref
    3147         try:
    3148             geo_reference = Geo_reference(NetCDFObject=fid)
    3149         except AttributeError, e:
    3150             geo_reference = Geo_reference() #Default georef object
    3151 
    3152         xllcorner = geo_reference.get_xllcorner()
    3153         yllcorner = geo_reference.get_yllcorner()
    3154         zone = geo_reference.get_zone()
    3155     else:
    3156         zone = origin[0]
    3157         xllcorner = origin[1]
    3158         yllcorner = origin[2]
    3159 
    3160 
    3161     #Get quantity and reduce if applicable
    3162     if verbose: print 'Reading quantity %s' %quantity
    3163 
    3164     if quantity.lower() == 'depth':
    3165         q = fid.variables['stage'][:] - fid.variables['elevation'][:]
    3166     else:
    3167         q = fid.variables[quantity][:]
    3168 
    3169 
    3170     if len(q.shape) == 2:
    3171         if verbose: print 'Reducing quantity %s' %quantity
    3172         q_reduced = zeros( number_of_points, Float )
    3173 
    3174         for k in range(number_of_points):
    3175             q_reduced[k] = reduction( q[:,k] )
    3176 
    3177         q = q_reduced
    3178 
    3179     #Now q has dimension: number_of_points
    3180 
    3181     #Create grid and update xll/yll corner and x,y
    3182     if verbose: print 'Creating grid'
    3183     ncols = int((xmax-xmin)/cellsize)+1
    3184     nrows = int((ymax-ymin)/cellsize)+1
    3185 
    3186     newxllcorner = xmin+xllcorner
    3187     newyllcorner = ymin+yllcorner
    3188 
    3189     x = x+xllcorner-newxllcorner
    3190     y = y+yllcorner-newyllcorner
    3191 
    3192     vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
    3193     assert len(vertex_points.shape) == 2
    3194 
    3195 
    3196     from Numeric import zeros, Float
    3197     grid_points = zeros ( (ncols*nrows, 2), Float )
    3198 
    3199 
    3200     for i in xrange(nrows):
    3201         yg = (nrows-i)*cellsize # this will flip the order of the y values
    3202         for j in xrange(ncols):
    3203             xg = j*cellsize
    3204             k = i*ncols + j
    3205 
    3206             grid_points[k,0] = xg
    3207             grid_points[k,1] = yg
    3208 
    3209     #Interpolate
    3210     from least_squares import Interpolation
    3211    
    3212 
    3213     #FIXME: This should be done with precrop = True (?), otherwise it'll
    3214     #take forever. With expand_search  set to False, some grid points might
    3215     #miss out....
    3216 
    3217     interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
    3218                            precrop = False, expand_search = True,
    3219                            verbose = verbose)
    3220 
    3221     #Interpolate using quantity values
    3222     if verbose: print 'Interpolating'
    3223     grid_values = interp.interpolate(q).flat
    3224     grid_values = reshape(grid_values,(nrows, ncols))
    3225 
    3226 
    3227     # setup header information
    3228     header['datum'] = '"' + datum + '"'
    3229     # FIXME The use of hardwired UTM and zone number needs to be made optional
    3230     # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
    3231     header['projection'] = '"UTM-' + str(zone) + '"'
    3232     header['coordinatetype'] = 'EN'
    3233     if header['coordinatetype'] == 'LL':
    3234         header['longitude'] = str(newxllcorner)
    3235         header['latitude'] = str(newyllcorner)
    3236     elif header['coordinatetype'] == 'EN':
    3237         header['eastings'] = str(newxllcorner)
    3238         header['northings'] = str(newyllcorner)
    3239     header['nullcellvalue'] = str(NODATA_value)
    3240     header['xdimension'] = str(cellsize)
    3241     header['ydimension'] = str(cellsize)
    3242     header['value'] = '"' + quantity + '"'
    3243 
    3244 
    3245     #Write
    3246     if verbose: print 'Writing %s' %ersfile
    3247     ermapper_grids.write_ermapper_grid(ersfile,grid_values, header)
    3248 
    3249     fid.close()   
    3250 
    3251 ##    ascid = open(ascfile, 'w')
    3252 ##
    3253 ##    ascid.write('ncols         %d\n' %ncols)
    3254 ##    ascid.write('nrows         %d\n' %nrows)
    3255 ##    ascid.write('xllcorner     %d\n' %newxllcorner)
    3256 ##    ascid.write('yllcorner     %d\n' %newyllcorner)
    3257 ##    ascid.write('cellsize      %f\n' %cellsize)
    3258 ##    ascid.write('NODATA_value  %d\n' %NODATA_value)
    3259 ##
    3260 ##
    3261 ##    #Get bounding polygon from mesh
    3262 ##    P = interp.mesh.get_boundary_polygon()
    3263 ##    inside_indices = inside_polygon(grid_points, P)
    3264 ##
    3265 ##    for i in range(nrows):
    3266 ##        if verbose and i%((nrows+10)/10)==0:
    3267 ##            print 'Doing row %d of %d' %(i, nrows)
    3268 ##
    3269 ##        for j in range(ncols):
    3270 ##            index = (nrows-i-1)*ncols+j
    3271 ##
    3272 ##            if sometrue(inside_indices == index):
    3273 ##                ascid.write('%f ' %grid_values[index])
    3274 ##            else:
    3275 ##                ascid.write('%d ' %NODATA_value)
    3276 ##
    3277 ##        ascid.write('\n')
    3278 ##
    3279 ##    #Close
    3280 ##    ascid.close()
    3281 ##    fid.close()
    3282 
Note: See TracChangeset for help on using the changeset viewer.