Changeset 2520
- Timestamp:
- Mar 10, 2006, 9:07:22 AM (19 years ago)
- Location:
- inundation/pyvolution
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/pyvolution/data_manager.py
r2518 r2520 1200 1200 totalnopoints = nrows*ncols 1201 1201 1202 # calculating number of NODATA_values for each row 1202 # calculating number of NODATA_values for each row in clipped region 1203 #FIXME: use array operations to do faster 1203 1204 nn = 0 1205 indices = [] 1204 1206 for i in range(nrows): 1205 1207 y = (nrows-i)*cellsize + yllcorner … … 1209 1211 northing_min <= y <= northing_max and \ 1210 1212 dem_elevation_r[i,j] == NODATA_value: 1211 nn += 1 1213 nn += 1 1214 1215 if easting_min <= x <= easting_max and \ 1216 northing_min <= y <= northing_max: 1217 indices.append([i,j]) 1218 1219 # finding indices of the clipped region 1220 if len(indices) > 0: 1221 firstpoint = indices[0] 1222 index1 = firstpoint[0] 1223 lastpoint = indices[len(indices)-1] 1224 index2 = lastpoint[0]+1 #? 1225 else: 1226 index1 = 1 1227 index2 = ncols 1212 1228 1213 1229 # dimension definitions … … 1222 1238 print 'There are %d values in the clipped elevation' %clippednopoints 1223 1239 print 'There are %d NODATA_values in the clipped elevation' %nn 1224 1240 1225 1241 outfile.createDimension('number_of_points', nopoints) 1226 1242 outfile.createDimension('number_of_dimensions', 2) #This is 2d data … … 1235 1251 elevation = outfile.variables['elevation'] 1236 1252 1253 1237 1254 #Store data 1238 1255 global_index = 0 … … 1242 1259 1243 1260 lower_index = global_index 1244 1245 v = dem_elevation_r[i, :]1261 1262 v = dem_elevation_r[i,index1:index2] 1246 1263 no_NODATA = sum(v == NODATA_value) 1247 1264 if no_NODATA > 0: … … 1249 1266 else: 1250 1267 newcols = ncols_in_bounding_box 1251 1268 1252 1269 telev = zeros(newcols, Float) 1253 1270 tpoints = zeros((newcols, 2), Float) -
inundation/pyvolution/test_data_manager.py
r2514 r2520 817 817 fid.close() 818 818 819 nn = ref_elevation.count(NODATA_value) 820 819 821 820 #Write prj file with metadata 822 821 metafilename = root+'.prj' … … 867 866 os.remove(root + '.prj') 868 867 868 def test_dem2pts_bounding_box_Nullvalues(self): 869 """Test conversion from dem in ascii format to native NetCDF xya format 870 """ 871 872 import time, os 873 from Numeric import array, zeros, allclose, Float, concatenate 874 from Scientific.IO.NetCDF import NetCDFFile 875 876 #Write test asc file 877 root = 'demtest' 878 879 filename = root+'.asc' 880 fid = open(filename, 'w') 881 fid.write("""ncols 5 882 nrows 6 883 xllcorner 2000.5 884 yllcorner 3000.5 885 cellsize 25 886 NODATA_value -9999 887 """) 888 #Create linear function 889 890 ref_points = [] 891 ref_elevation = [] 892 new_ref_pts1 = [] 893 new_ref_elev1 = [] 894 NODATA_value = -9999 895 for i in range(6): 896 y = (6-i)*25.0 897 for j in range(5): 898 x = j*25.0 899 z = x+2*y 900 if j == 4: z = NODATA_value # column 901 if i == 2 and j == 2: z = NODATA_value # random 902 if i == 5 and j == 1: z = NODATA_value 903 if i == 1: z = NODATA_value # row 904 if i == 3 and j == 1: z = NODATA_value # two pts/row 905 if i == 3 and j == 3: z = NODATA_value 906 907 if z <> NODATA_value: 908 new_ref_elev1.append(z) 909 new_ref_pts1.append( [x,y] ) 910 911 ref_points.append( [x,y] ) 912 ref_elevation.append(z) 913 fid.write('%f ' %z) 914 fid.write('\n') 915 916 fid.close() 917 918 #Write prj file with metadata 919 metafilename = root+'.prj' 920 fid = open(metafilename, 'w') 921 922 923 fid.write("""Projection UTM 924 Zone 56 925 Datum WGS84 926 Zunits NO 927 Units METERS 928 Spheroid WGS84 929 Xshift 0.0000000000 930 Yshift 10000000.0000000000 931 Parameters 932 """) 933 fid.close() 934 935 #Convert to NetCDF pts 936 convert_dem_from_ascii2netcdf(root) 937 dem2pts(root, easting_min=2010.0, easting_max=2110.0, 938 northing_min=3035.0, northing_max=3125.5) 939 940 #Check contents 941 #Get NetCDF 942 fid = NetCDFFile(root+'.pts', 'r') 943 944 # Get the variables 945 #print fid.variables.keys() 946 points = fid.variables['points'] 947 elevation = fid.variables['elevation'] 948 949 #Check values 950 assert fid.xllcorner[0] == 2010.0 951 assert fid.yllcorner[0] == 3035.0 952 953 #create new reference points 954 ref_points = [] 955 ref_elevation = [] 956 new_ref_pts2 = [] 957 new_ref_elev2 = [] 958 for i in range(4): 959 y = (4-i)*25.0 + 25.0 960 y_new = y + 3000.5 - 3035.0 961 for j in range(4): 962 x = j*25.0 + 25.0 963 x_new = x + 2000.5 - 2010.0 964 z = x+2*y 965 966 if j == 3: z = NODATA_value # column 967 if i == 1 and j == 1: z = NODATA_value # random 968 if i == 4 and j == 0: z = NODATA_value 969 if i == 0: z = NODATA_value # row 970 if i == 2 and j == 0: z = NODATA_value # two pts/row 971 if i == 2 and j == 2: z = NODATA_value 972 973 if z <> NODATA_value: 974 new_ref_elev2.append(z) 975 new_ref_pts2.append( [x_new,y_new] ) 976 977 978 ref_points.append( [x_new,y_new] ) 979 ref_elevation.append(z) 980 981 #print points[:] 982 #print ref_points 983 #assert allclose(points, ref_points) 984 985 #print attributes[:] 986 #print ref_elevation 987 #assert allclose(elevation, ref_elevation) 988 989 990 assert len(points) == len(new_ref_pts2), 'length of returned points not correct' 991 assert allclose(points, new_ref_pts2), 'points do not align' 992 993 #print 'elevation', elevation[:] 994 assert len(elevation) == len(new_ref_elev2), 'length of returned elevation not correct' 995 assert allclose(elevation, new_ref_elev2), 'elevations do not align' 996 #Cleanup 997 fid.close() 998 999 1000 os.remove(root + '.pts') 1001 os.remove(root + '.dem') 1002 os.remove(root + '.asc') 1003 os.remove(root + '.prj') 869 1004 870 1005
Note: See TracChangeset
for help on using the changeset viewer.