Changeset 1992
- Timestamp:
- Nov 2, 2005, 10:20:13 AM (18 years ago)
- Location:
- inundation/pyvolution
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/pyvolution/data_manager.py
r1985 r1992 49 49 """ 50 50 51 from Numeric import concatenate 51 from Numeric import concatenate, array, Float 52 52 53 53 from coordinate_transforms.geo_reference import Geo_reference … … 2839 2839 outfile.close() 2840 2840 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 asc_csiro2sww(bath_file_in, verbose=False): 3059 pass 3060 3061 def _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 6 6 from Numeric import zeros, array, allclose, Float 7 7 from util import mean 8 import tempfile 8 9 9 10 from data_manager import * … … 2653 2654 verbose = True) 2654 2655 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 2669 nrows 4 2670 xllcorner 2000.5 2671 yllcorner 3000.5 2672 cellsize 25 2673 NODATA_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 2657 2689 2658 2690 #-------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.