Changeset 1985
 Timestamp:
 Oct 31, 2005, 10:38:30 PM (18 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

inundation/pyvolution/data_manager.py
r1984 r1985 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 31212855 nrows 18002856 xllcorner 7220002857 yllcorner 58930002858 cellsize 252859 NODATA_value 99992860 138.3698 137.4194 136.5062 135.5558 ..........2861 2862 Also write accompanying file with same basename_in but extension .prj2863 used to fix the UTM zone, datum, false northings and eastings.2864 2865 The prj format is assumed to be as2866 2867 Projection UTM2868 Zone 562869 Datum WGS842870 Zunits NO2871 Units METERS2872 Spheroid WGS842873 Xshift 0.00000000002874 Yshift 10000000.00000000002875 Parameters2876 2877 2878 if quantity is given, out values from quantity otherwise default to2879 elevation2880 2881 if timestep (an index) is given, output quantity at that timestep2882 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 sometrue2888 from utilities.polygon import inside_polygon2889 2890 #FIXME: Should be variable2891 datum = 'WGS84'2892 false_easting = 5000002893 false_northing = 100000002894 2895 if quantity is None:2896 quantity = 'elevation'2897 2898 if reduction is None:2899 reduction = max2900 2901 if basename_out is None:2902 basename_out = basename_in + '_%s' %quantity2903 2904 swwfile = basename_in + '.sww'2905 ascfile = basename_out + '.asc'2906 prjfile = basename_out + '.prj'2907 2908 2909 if verbose: print 'Reading from %s' %swwfile2910 #Read sww file2911 from Scientific.IO.NetCDF import NetCDFFile2912 fid = NetCDFFile(swwfile)2913 2914 #Get extent and reference2915 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_reference2927 #sww files don't have to have a geo_ref2928 try:2929 geo_reference = Geo_reference(NetCDFObject=fid)2930 except AttributeError, e:2931 geo_reference = Geo_reference() #Default georef object2932 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 applicable2943 if verbose: print 'Reading quantity %s' %quantity2944 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' %quantity2953 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_reduced2959 2960 #Now q has dimension: number_of_points2961 2962 #Create grid and update xll/yll corner and x,y2963 if verbose: print 'Creating grid'2964 ncols = int((xmaxxmin)/cellsize)+12965 nrows = int((ymaxymin)/cellsize)+12966 2967 newxllcorner = xmin+xllcorner2968 newyllcorner = ymin+yllcorner2969 2970 x = x+xllcornernewxllcorner2971 y = y+yllcornernewyllcorner2972 2973 vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)2974 assert len(vertex_points.shape) == 22975 2976 2977 from Numeric import zeros, Float2978 grid_points = zeros ( (ncols*nrows, 2), Float )2979 2980 2981 for i in xrange(nrows):2982 yg = i*cellsize2983 for j in xrange(ncols):2984 xg = j*cellsize2985 k = i*ncols + j2986 2987 grid_points[k,0] = xg2988 grid_points[k,1] = yg2989 2990 #Interpolate2991 from least_squares import Interpolation2992 2993 2994 #FIXME: This should be done with precrop = True, otherwise it'll2995 #take forever. With expand_search set to False, some grid points might2996 #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 values3003 if verbose: print 'Interpolating'3004 grid_values = interp.interpolate(q).flat3005 3006 #Write3007 #Write prj file3008 if verbose: print 'Writing %s' %prjfile3009 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' %ascfile3024 NODATA_value = 99993025 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 mesh3037 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 = (nrowsi1)*ncols+j3046 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 #Close3055 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 31213072 nrows 18003073 xllcorner 7220003074 yllcorner 58930003075 cellsize 253076 NODATA_value 99993077 138.3698 137.4194 136.5062 135.5558 ..........3078 3079 Also write accompanying file with same basename_in but extension .prj3080 used to fix the UTM zone, datum, false northings and eastings.3081 3082 The prj format is assumed to be as3083 3084 Projection UTM3085 Zone 563086 Datum WGS843087 Zunits NO3088 Units METERS3089 Spheroid WGS843090 Xshift 0.00000000003091 Yshift 10000000.00000000003092 Parameters3093 3094 3095 if quantity is given, out values from quantity otherwise default to3096 elevation3097 3098 if timestep (an index) is given, output quantity at that timestep3099 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, sometrue3104 import ermapper_grids3105 from utilities.polygon import inside_polygon3106 3107 header = {}3108 false_easting = 5000003109 false_northing = 100000003110 NODATA_value = 03111 3112 if quantity is None:3113 quantity = 'elevation'3114 3115 if reduction is None:3116 reduction = max3117 3118 if basename_out is None:3119 basename_out = basename_in + '_%s' %quantity3120 3121 swwfile = basename_in + '.sww'3122 # Note the use of a .ers extension is optional (write_ermapper_grid will3123 # deal with either option3124 ersfile = basename_out3125 ## prjfile = basename_out + '.prj'3126 3127 3128 if verbose: print 'Reading from %s' %swwfile3129 #Read sww file3130 from Scientific.IO.NetCDF import NetCDFFile3131 fid = NetCDFFile(swwfile)3132 3133 #Get extent and reference3134 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_reference3146 #sww files don't have to have a geo_ref3147 try:3148 geo_reference = Geo_reference(NetCDFObject=fid)3149 except AttributeError, e:3150 geo_reference = Geo_reference() #Default georef object3151 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 applicable3162 if verbose: print 'Reading quantity %s' %quantity3163 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' %quantity3172 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_reduced3178 3179 #Now q has dimension: number_of_points3180 3181 #Create grid and update xll/yll corner and x,y3182 if verbose: print 'Creating grid'3183 ncols = int((xmaxxmin)/cellsize)+13184 nrows = int((ymaxymin)/cellsize)+13185 3186 newxllcorner = xmin+xllcorner3187 newyllcorner = ymin+yllcorner3188 3189 x = x+xllcornernewxllcorner3190 y = y+yllcornernewyllcorner3191 3192 vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)3193 assert len(vertex_points.shape) == 23194 3195 3196 from Numeric import zeros, Float3197 grid_points = zeros ( (ncols*nrows, 2), Float )3198 3199 3200 for i in xrange(nrows):3201 yg = (nrowsi)*cellsize # this will flip the order of the y values3202 for j in xrange(ncols):3203 xg = j*cellsize3204 k = i*ncols + j3205 3206 grid_points[k,0] = xg3207 grid_points[k,1] = yg3208 3209 #Interpolate3210 from least_squares import Interpolation3211 3212 3213 #FIXME: This should be done with precrop = True (?), otherwise it'll3214 #take forever. With expand_search set to False, some grid points might3215 #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 values3222 if verbose: print 'Interpolating'3223 grid_values = interp.interpolate(q).flat3224 grid_values = reshape(grid_values,(nrows, ncols))3225 3226 3227 # setup header information3228 header['datum'] = '"' + datum + '"'3229 # FIXME The use of hardwired UTM and zone number needs to be made optional3230 # 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 #Write3246 if verbose: print 'Writing %s' %ersfile3247 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 mesh3262 ## 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 = (nrowsi1)*ncols+j3271 ##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 ## #Close3280 ## ascid.close()3281 ## fid.close()3282
Note: See TracChangeset
for help on using the changeset viewer.