Changeset 2553
 Timestamp:
 Mar 16, 2006, 5:35:18 PM (18 years ago)
 Location:
 inundation/pyvolution
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

inundation/pyvolution/data_manager.py
r2541 r2553 1206 1206 i1_0 = 0 1207 1207 j1_0 = 0 1208 thisj = 0 1209 thisi = 0 1208 1210 for i in range(nrows): 1209 y = (nrowsi )*cellsize + yllcorner1211 y = (nrowsi1)*cellsize + yllcorner 1210 1212 for j in range(ncols): 1211 1213 x = j*cellsize + xllcorner 1212 1214 if easting_min <= x <= easting_max and \ 1213 1215 northing_min <= y <= northing_max: 1216 thisj = j 1217 thisi = i 1214 1218 if dem_elevation_r[i,j] == NODATA_value: nn += 1 1215 1219 1216 1220 if k == 0: 1217 1221 i1_0 = i 1218 j1_0 = j 1222 j1_0 = j 1219 1223 k += 1 1220 1224 1221 index1 = j1_0 1222 index2 = j 1223 index3 = i1_0 1224 index4 = i 1225 1226 1227 index2 += 1 1228 nrows2 = index4  index3 1229 ncols2 = index2  index1 + 1 1225 index1 = j1_0 1226 index2 = thisj 1230 1227 1231 1228 # dimension definitions … … 1233 1230 ncols_in_bounding_box = int(round((easting_maxeasting_min)/cellsize)) 1234 1231 1235 clippednopoints = nrows_in_bounding_box*ncols_in_bounding_box1232 clippednopoints = (thisi+1i1_0)*(thisj+1j1_0) 1236 1233 nopoints = clippednopointsnn 1237 1234 1235 clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1] 1236 1238 1237 if verbose and nn > 0: 1239 1238 print 'There are %d values in the elevation' %totalnopoints 1240 #print 'There are %d values in the clipped elevation' %clippednopoints1239 print 'There are %d values in the clipped elevation' %clippednopoints 1241 1240 print 'There are %d NODATA_values in the clipped elevation' %nn 1242 1241 … … 1253 1252 elevation = outfile.variables['elevation'] 1254 1253 1255 1254 lenv = index2index1+1 1256 1255 #Store data 1257 1256 global_index = 0 1258 for i in range(nrows): 1257 #for i in range(nrows): 1258 for i in range(i1_0,thisi+1,1): 1259 1259 if verbose and i%((nrows+10)/10)==0: 1260 1260 print 'Processing row %d of %d' %(i, nrows) … … 1262 1262 lower_index = global_index 1263 1263 1264 v = dem_elevation_r[i,index1:index2 ]1264 v = dem_elevation_r[i,index1:index2+1] 1265 1265 no_NODATA = sum(v == NODATA_value) 1266 1266 if no_NODATA > 0: 1267 newcols = len (v) no_NODATA #ncols_in_bounding_box  no_NODATA1267 newcols = lenv  no_NODATA #ncols_in_bounding_box  no_NODATA 1268 1268 else: 1269 newcols = len(v) #ncols_in_bounding_box 1270 1271 #print 'here', len(v), no_NODATA, newcols, v 1272 #print 'here again', dem_elevation_r[i,:] 1269 newcols = lenv #ncols_in_bounding_box 1270 1273 1271 telev = zeros(newcols, Float) 1274 1272 tpoints = zeros((newcols, 2), Float) … … 1276 1274 local_index = 0 1277 1275 1278 y = (nrowsi)*cellsize + yllcorner 1279 for j in range(ncols): 1276 y = (nrowsi1)*cellsize + yllcorner 1277 #for j in range(ncols): 1278 for j in range(j1_0,index2+1,1): 1280 1279 1281 1280 x = j*cellsize + xllcorner 
inundation/pyvolution/test_data_manager.py
r2526 r2553 571 571 572 572 573 574 575 def test_dem2pts(self): 573 def test_dem2pts_bounding_box_v2(self): 576 574 """Test conversion from dem in ascii format to native NetCDF xya format 577 575 """ 578 576 579 577 import time, os 580 from Numeric import array, zeros, allclose, Float, concatenate 578 from Numeric import array, zeros, allclose, Float, concatenate, ones 581 579 from Scientific.IO.NetCDF import NetCDFFile 582 580 … … 586 584 filename = root+'.asc' 587 585 fid = open(filename, 'w') 588 fid.write("""ncols 5589 nrows 6590 xllcorner 2000 .5591 yllcorner 3000 .5592 cellsize 25586 fid.write("""ncols 10 587 nrows 10 588 xllcorner 2000 589 yllcorner 3000 590 cellsize 1 593 591 NODATA_value 9999 594 592 """) 595 #Create linear function 596 593 #Create linear function 597 594 ref_points = [] 598 595 ref_elevation = [] 599 for i in range(6): 600 y = (6i)*25.0 601 for j in range(5): 602 x = j*25.0 603 z = x+2*y 604 605 ref_points.append( [x,y] ) 596 x0 = 2000 597 y = 3010 598 yvec = range(10) 599 xvec = range(10) 600 z = 1 601 for i in range(10): 602 y = y  1 603 for j in range(10): 604 x = x0 + xvec[j] 605 z += 1 606 ref_points.append ([x,y]) 606 607 ref_elevation.append(z) 607 608 fid.write('%f ' %z) … … 609 610 610 611 fid.close() 612 613 #print 'sending pts', ref_points 614 #print 'sending elev', ref_elevation 611 615 612 616 #Write prj file with metadata … … 629 633 #Convert to NetCDF pts 630 634 convert_dem_from_ascii2netcdf(root) 631 dem2pts(root) 635 dem2pts(root, easting_min=2002.0, easting_max=2007.0, 636 northing_min=3003.0, northing_max=3006.0) 632 637 633 638 #Check contents … … 641 646 642 647 #Check values 643 644 #print points[:] 645 #print ref_points 646 assert allclose(points, ref_points) 648 assert fid.xllcorner[0] == 2002.0 649 assert fid.yllcorner[0] == 3003.0 650 651 #create new reference points 652 newz = [] 653 newz[0:5] = ref_elevation[32]*ones(6) + range(6) 654 newz[6:11] = ref_elevation[42]*ones(6) + range(6) 655 newz[12:17] = ref_elevation[52]*ones(6) + range(6) 656 newz[18:23] = ref_elevation[62]*ones(6) + range(6) 657 ref_elevation = [] 658 ref_elevation = newz 659 ref_points = [] 660 new_ref_points = [] 661 x0 = 2002 662 y = 3007 663 yvec = range(4) 664 xvec = range(6) 665 for i in range(4): 666 y = y  1 667 ynew = y  3003.0 668 for j in range(6): 669 x = x0 + xvec[j] 670 xnew = x  2002.0 671 ref_points.append([x,y]) 672 new_ref_points.append ([xnew,ynew]) 673 674 #print 'bounding box' 675 #print 'new ref pts', new_ref_points[:] 676 #print 'received pts', points[:] 677 678 assert allclose(points, new_ref_points) 647 679 648 680 #print attributes[:] 649 #print ref_elevation 681 #print 'ref_elev', ref_elevation[:] 682 #print 'received elev', elevation[:] 650 683 assert allclose(elevation, ref_elevation) 651 684 … … 660 693 661 694 662 663 def test_dem2pts_bounding_box(self): 695 def test_dem2pts_bounding_box_removeNullvalues_v2(self): 664 696 """Test conversion from dem in ascii format to native NetCDF xya format 665 697 """ 666 698 667 699 import time, os 668 from Numeric import array, zeros, allclose, Float, concatenate 700 from Numeric import array, zeros, allclose, Float, concatenate, ones 669 701 from Scientific.IO.NetCDF import NetCDFFile 670 702 … … 674 706 filename = root+'.asc' 675 707 fid = open(filename, 'w') 676 fid.write("""ncols 5677 nrows 6678 xllcorner 2000 .5679 yllcorner 3000 .5680 cellsize 25708 fid.write("""ncols 10 709 nrows 10 710 xllcorner 2000 711 yllcorner 3000 712 cellsize 1 681 713 NODATA_value 9999 682 714 """) 683 #Create linear function 684 715 #Create linear function 685 716 ref_points = [] 686 717 ref_elevation = [] 687 for i in range(6): 688 y = (6i)*25.0 689 for j in range(5): 690 x = j*25.0 691 z = x+2*y 692 693 ref_points.append( [x,y] ) 694 ref_elevation.append(z) 695 fid.write('%f ' %z) 718 x0 = 2000 719 y = 3010 720 yvec = range(10) 721 xvec = range(10) 722 z = range(100) 723 NODATA_value = 9999 724 count = 1 725 for i in range(10): 726 y = y  1 727 for j in range(10): 728 x = x0 + xvec[j] 729 ref_points.append ([x,y]) 730 count += 1 731 if j == 4: z[count] = NODATA_value #row 732 if j == 8: z[count] = NODATA_value #row 733 if i == 9: z[count] = NODATA_value #column 734 ref_elevation.append( z[count] ) 735 fid.write('%f ' %z[count]) 696 736 fid.write('\n') 697 737 698 fid.close() 738 fid.close() 739 740 #print 'sending elev', ref_elevation 699 741 700 742 #Write prj file with metadata … … 717 759 #Convert to NetCDF pts 718 760 convert_dem_from_ascii2netcdf(root) 719 dem2pts(root, easting_min=20 10.0, easting_max=2110.0,720 northing_min=30 35.0, northing_max=3125.5)761 dem2pts(root, easting_min=2002.0, easting_max=2007.0, 762 northing_min=3003.0, northing_max=3006.0) 721 763 722 764 #Check contents … … 730 772 731 773 #Check values 732 assert fid.xllcorner[0] == 20 10.0733 assert fid.yllcorner[0] == 30 35.0774 assert fid.xllcorner[0] == 2002.0 775 assert fid.yllcorner[0] == 3003.0 734 776 735 777 #create new reference points 778 newz = [] 779 newz[0:1] = ref_elevation[32]*ones(2) + range(2) 780 newz[2:4] = ref_elevation[32]*ones(3) + range(3,6) 781 newz[5:6] = ref_elevation[42]*ones(2) + range(2) 782 newz[7:9] = ref_elevation[42]*ones(3) + range(3,6) 783 newz[10:11] = ref_elevation[52]*ones(2) + range(2) 784 newz[12:14] = ref_elevation[52]*ones(3) + range(3,6) 785 newz[15:16] = ref_elevation[62]*ones(2) + range(2) 786 newz[17:19] = ref_elevation[62]*ones(3) + range(3,6) 787 788 ref_elevation = [] 789 ref_elevation = newz 736 790 ref_points = [] 737 ref_elevation = [] 791 new_ref_points = [] 792 x0 = 2002 793 y = 3007 794 yvec = range(4) 795 xvec = range(6) 738 796 for i in range(4): 739 y = (4i)*25.0 + 25.0 740 y_new = y + 3000.5  3035.0 741 for j in range(4): 742 x = j*25.0 + 25.0 743 x_new = x + 2000.5  2010.0 744 z = x+2*y 745 746 ref_points.append( [x_new,y_new] ) 747 ref_elevation.append(z) 748 749 #print points[:] 750 #print ref_points 751 assert allclose(points, ref_points) 752 753 #print attributes[:] 754 #print ref_elevation 797 y = y  1 798 ynew = y  3003.0 799 for j in range(6): 800 x = x0 + xvec[j] 801 xnew = x  2002.0 802 if j <> 2: 803 ref_points.append([x,y]) 804 new_ref_points.append ([xnew,ynew]) 805 806 assert allclose(points, new_ref_points) 755 807 assert allclose(elevation, ref_elevation) 756 808 … … 763 815 os.remove(root + '.asc') 764 816 os.remove(root + '.prj') 765 766 def test_dem2pts_remove_Nullvalues(self):767 """Test conversion from dem in ascii format to native NetCDF xya format768 """769 770 import time, os771 from Numeric import array, zeros, allclose, Float, concatenate772 from Scientific.IO.NetCDF import NetCDFFile773 774 #Write test asc file775 root = 'demtest'776 777 filename = root+'.asc'778 fid = open(filename, 'w')779 fid.write("""ncols 5780 nrows 6781 xllcorner 2000.5782 yllcorner 3000.5783 cellsize 25784 NODATA_value 9999785 """)786 #Create linear function787 # ref_ will write all the values788 # new_ref_ will write the values except for NODATA_values789 ref_points = []790 ref_elevation = []791 new_ref_pts = []792 new_ref_elev = []793 NODATA_value = 9999794 for i in range(6):795 y = (6i)*25.0796 for j in range(5):797 x = j*25.0798 z = x+2*y799 if j == 4: z = NODATA_value # column800 if i == 2 and j == 2: z = NODATA_value # random801 if i == 5 and j == 1: z = NODATA_value802 if i == 1: z = NODATA_value # row803 if i == 3 and j == 1: z = NODATA_value # two pts/row804 if i == 3 and j == 3: z = NODATA_value805 806 807 if z <> NODATA_value:808 new_ref_elev.append(z)809 new_ref_pts.append( [x,y] )810 811 ref_points.append( [x,y] )812 ref_elevation.append(z)813 814 fid.write('%f ' %z)815 fid.write('\n')816 817 fid.close()818 819 820 #Write prj file with metadata821 metafilename = root+'.prj'822 fid = open(metafilename, 'w')823 824 825 fid.write("""Projection UTM826 Zone 56827 Datum WGS84828 Zunits NO829 Units METERS830 Spheroid WGS84831 Xshift 0.0000000000832 Yshift 10000000.0000000000833 Parameters834 """)835 fid.close()836 837 #Convert to NetCDF pts838 convert_dem_from_ascii2netcdf(root)839 dem2pts(root)840 841 #Check contents842 #Get NetCDF843 fid = NetCDFFile(root+'.pts', 'r')844 845 # Get the variables846 #print fid.variables.keys()847 points = fid.variables['points']848 elevation = fid.variables['elevation']849 850 #Check values851 #print 'points', points[:]852 assert len(points) == len(new_ref_pts), 'length of returned points not correct'853 assert allclose(points, new_ref_pts), 'points do not align'854 855 #print 'elevation', elevation[:]856 assert len(elevation) == len(new_ref_elev), 'length of returned elevation not correct'857 assert allclose(elevation, new_ref_elev), 'elevations do not align'858 859 #Cleanup860 fid.close()861 862 863 os.remove(root + '.pts')864 os.remove(root + '.dem')865 os.remove(root + '.asc')866 os.remove(root + '.prj')867 868 def test_dem2pts_bounding_box_Nullvalues(self):869 """Test conversion from dem in ascii format to native NetCDF xya format870 """871 872 import time, os873 from Numeric import array, zeros, allclose, Float, concatenate874 from Scientific.IO.NetCDF import NetCDFFile875 876 #Write test asc file877 root = 'demtest'878 879 filename = root+'.asc'880 fid = open(filename, 'w')881 fid.write("""ncols 5882 nrows 6883 xllcorner 2000.5884 yllcorner 3000.5885 cellsize 25886 NODATA_value 9999887 """)888 #Create linear function889 890 ref_points = []891 ref_elevation = []892 new_ref_pts1 = []893 new_ref_elev1 = []894 NODATA_value = 9999895 for i in range(6):896 y = (6i)*25.0897 for j in range(5):898 x = j*25.0899 z = x+2*y900 if j == 4: z = NODATA_value # column901 if i == 2 and j == 2: z = NODATA_value # random902 if i == 5 and j == 1: z = NODATA_value903 if i == 1: z = NODATA_value # row904 if i == 3 and j == 1: z = NODATA_value # two pts/row905 if i == 3 and j == 3: z = NODATA_value906 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 metadata919 metafilename = root+'.prj'920 fid = open(metafilename, 'w')921 922 923 fid.write("""Projection UTM924 Zone 56925 Datum WGS84926 Zunits NO927 Units METERS928 Spheroid WGS84929 Xshift 0.0000000000930 Yshift 10000000.0000000000931 Parameters932 """)933 fid.close()934 935 #Convert to NetCDF pts936 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 contents941 #Get NetCDF942 fid = NetCDFFile(root+'.pts', 'r')943 944 # Get the variables945 #print fid.variables.keys()946 points = fid.variables['points']947 elevation = fid.variables['elevation']948 949 #Check values950 assert fid.xllcorner[0] == 2010.0951 assert fid.yllcorner[0] == 3035.0952 953 #create new reference points954 ref_points = []955 ref_elevation = []956 new_ref_pts2 = []957 new_ref_elev2 = []958 for i in range(4):959 y = (4i)*25.0 + 25.0960 y_new = y + 3000.5  3035.0961 for j in range(4):962 x = j*25.0 + 25.0963 x_new = x + 2000.5  2010.0964 z = x+2*y965 966 if j == 3: z = NODATA_value # column967 if i == 1 and j == 1: z = NODATA_value # random968 if i == 4 and j == 0: z = NODATA_value969 if i == 0: z = NODATA_value # row970 if i == 2 and j == 0: z = NODATA_value # two pts/row971 if i == 2 and j == 2: z = NODATA_value972 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_points983 #assert allclose(points, ref_points)984 985 #print attributes[:]986 #print ref_elevation987 #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 #Cleanup997 fid.close()998 999 1000 os.remove(root + '.pts')1001 os.remove(root + '.dem')1002 os.remove(root + '.asc')1003 os.remove(root + '.prj')1004 1005 1006 817 1007 818 … … 4041 3852 runner.run(suite) 4042 3853 3854 3855 3856 3857 def test_dem2pts(self): 3858 """Test conversion from dem in ascii format to native NetCDF xya format 3859 """ 3860 3861 import time, os 3862 from Numeric import array, zeros, allclose, Float, concatenate 3863 from Scientific.IO.NetCDF import NetCDFFile 3864 3865 #Write test asc file 3866 root = 'demtest' 3867 3868 filename = root+'.asc' 3869 fid = open(filename, 'w') 3870 fid.write("""ncols 5 3871 nrows 6 3872 xllcorner 2000.5 3873 yllcorner 3000.5 3874 cellsize 25 3875 NODATA_value 9999 3876 """) 3877 #Create linear function 3878 3879 ref_points = [] 3880 ref_elevation = [] 3881 for i in range(6): 3882 y = (6i)*25.0 3883 for j in range(5): 3884 x = j*25.0 3885 z = x+2*y 3886 3887 ref_points.append( [x,y] ) 3888 ref_elevation.append(z) 3889 fid.write('%f ' %z) 3890 fid.write('\n') 3891 3892 fid.close() 3893 3894 #Write prj file with metadata 3895 metafilename = root+'.prj' 3896 fid = open(metafilename, 'w') 3897 3898 3899 fid.write("""Projection UTM 3900 Zone 56 3901 Datum WGS84 3902 Zunits NO 3903 Units METERS 3904 Spheroid WGS84 3905 Xshift 0.0000000000 3906 Yshift 10000000.0000000000 3907 Parameters 3908 """) 3909 fid.close() 3910 3911 #Convert to NetCDF pts 3912 convert_dem_from_ascii2netcdf(root) 3913 dem2pts(root) 3914 3915 #Check contents 3916 #Get NetCDF 3917 fid = NetCDFFile(root+'.pts', 'r') 3918 3919 # Get the variables 3920 #print fid.variables.keys() 3921 points = fid.variables['points'] 3922 elevation = fid.variables['elevation'] 3923 3924 #Check values 3925 3926 #print points[:] 3927 #print ref_points 3928 assert allclose(points, ref_points) 3929 3930 #print attributes[:] 3931 #print ref_elevation 3932 assert allclose(elevation, ref_elevation) 3933 3934 #Cleanup 3935 fid.close() 3936 3937 3938 os.remove(root + '.pts') 3939 os.remove(root + '.dem') 3940 os.remove(root + '.asc') 3941 os.remove(root + '.prj') 3942 3943 3944 3945 def test_dem2pts_bounding_box(self): 3946 """Test conversion from dem in ascii format to native NetCDF xya format 3947 """ 3948 3949 import time, os 3950 from Numeric import array, zeros, allclose, Float, concatenate 3951 from Scientific.IO.NetCDF import NetCDFFile 3952 3953 #Write test asc file 3954 root = 'demtest' 3955 3956 filename = root+'.asc' 3957 fid = open(filename, 'w') 3958 fid.write("""ncols 5 3959 nrows 6 3960 xllcorner 2000.5 3961 yllcorner 3000.5 3962 cellsize 25 3963 NODATA_value 9999 3964 """) 3965 #Create linear function 3966 3967 ref_points = [] 3968 ref_elevation = [] 3969 for i in range(6): 3970 y = (6i)*25.0 3971 for j in range(5): 3972 x = j*25.0 3973 z = x+2*y 3974 3975 ref_points.append( [x,y] ) 3976 ref_elevation.append(z) 3977 fid.write('%f ' %z) 3978 fid.write('\n') 3979 3980 fid.close() 3981 3982 #Write prj file with metadata 3983 metafilename = root+'.prj' 3984 fid = open(metafilename, 'w') 3985 3986 3987 fid.write("""Projection UTM 3988 Zone 56 3989 Datum WGS84 3990 Zunits NO 3991 Units METERS 3992 Spheroid WGS84 3993 Xshift 0.0000000000 3994 Yshift 10000000.0000000000 3995 Parameters 3996 """) 3997 fid.close() 3998 3999 #Convert to NetCDF pts 4000 convert_dem_from_ascii2netcdf(root) 4001 dem2pts(root, easting_min=2010.0, easting_max=2110.0, 4002 northing_min=3035.0, northing_max=3125.5) 4003 4004 #Check contents 4005 #Get NetCDF 4006 fid = NetCDFFile(root+'.pts', 'r') 4007 4008 # Get the variables 4009 #print fid.variables.keys() 4010 points = fid.variables['points'] 4011 elevation = fid.variables['elevation'] 4012 4013 #Check values 4014 assert fid.xllcorner[0] == 2010.0 4015 assert fid.yllcorner[0] == 3035.0 4016 4017 #create new reference points 4018 ref_points = [] 4019 ref_elevation = [] 4020 for i in range(4): 4021 y = (4i)*25.0 + 25.0 4022 y_new = y + 3000.5  3035.0 4023 for j in range(4): 4024 x = j*25.0 + 25.0 4025 x_new = x + 2000.5  2010.0 4026 z = x+2*y 4027 4028 ref_points.append( [x_new,y_new] ) 4029 ref_elevation.append(z) 4030 4031 #print points[:] 4032 #print ref_points 4033 assert allclose(points, ref_points) 4034 4035 #print attributes[:] 4036 #print ref_elevation 4037 assert allclose(elevation, ref_elevation) 4038 4039 #Cleanup 4040 fid.close() 4041 4042 4043 os.remove(root + '.pts') 4044 os.remove(root + '.dem') 4045 os.remove(root + '.asc') 4046 os.remove(root + '.prj') 4047 4048 4049 4050 def test_dem2pts_remove_Nullvalues(self): 4051 """Test conversion from dem in ascii format to native NetCDF xya format 4052 """ 4053 4054 import time, os 4055 from Numeric import array, zeros, allclose, Float, concatenate 4056 from Scientific.IO.NetCDF import NetCDFFile 4057 4058 #Write test asc file 4059 root = 'demtest' 4060 4061 filename = root+'.asc' 4062 fid = open(filename, 'w') 4063 fid.write("""ncols 5 4064 nrows 6 4065 xllcorner 2000.5 4066 yllcorner 3000.5 4067 cellsize 25 4068 NODATA_value 9999 4069 """) 4070 #Create linear function 4071 # ref_ will write all the values 4072 # new_ref_ will write the values except for NODATA_values 4073 ref_points = [] 4074 ref_elevation = [] 4075 new_ref_pts = [] 4076 new_ref_elev = [] 4077 NODATA_value = 9999 4078 for i in range(6): 4079 y = (6i)*25.0 4080 for j in range(5): 4081 x = j*25.0 4082 z = x+2*y 4083 if j == 4: z = NODATA_value # column 4084 if i == 2 and j == 2: z = NODATA_value # random 4085 if i == 5 and j == 1: z = NODATA_value 4086 if i == 1: z = NODATA_value # row 4087 if i == 3 and j == 1: z = NODATA_value # two pts/row 4088 if i == 3 and j == 3: z = NODATA_value 4089 4090 4091 if z <> NODATA_value: 4092 new_ref_elev.append(z) 4093 new_ref_pts.append( [x,y] ) 4094 4095 ref_points.append( [x,y] ) 4096 ref_elevation.append(z) 4097 4098 fid.write('%f ' %z) 4099 fid.write('\n') 4100 4101 fid.close() 4102 4103 4104 #Write prj file with metadata 4105 metafilename = root+'.prj' 4106 fid = open(metafilename, 'w') 4107 4108 4109 fid.write("""Projection UTM 4110 Zone 56 4111 Datum WGS84 4112 Zunits NO 4113 Units METERS 4114 Spheroid WGS84 4115 Xshift 0.0000000000 4116 Yshift 10000000.0000000000 4117 Parameters 4118 """) 4119 fid.close() 4120 4121 #Convert to NetCDF pts 4122 convert_dem_from_ascii2netcdf(root) 4123 dem2pts(root) 4124 4125 #Check contents 4126 #Get NetCDF 4127 fid = NetCDFFile(root+'.pts', 'r') 4128 4129 # Get the variables 4130 #print fid.variables.keys() 4131 points = fid.variables['points'] 4132 elevation = fid.variables['elevation'] 4133 4134 #Check values 4135 #print 'points', points[:] 4136 assert len(points) == len(new_ref_pts), 'length of returned points not correct' 4137 assert allclose(points, new_ref_pts), 'points do not align' 4138 4139 #print 'elevation', elevation[:] 4140 assert len(elevation) == len(new_ref_elev), 'length of returned elevation not correct' 4141 assert allclose(elevation, new_ref_elev), 'elevations do not align' 4142 4143 #Cleanup 4144 fid.close() 4145 4146 4147 os.remove(root + '.pts') 4148 os.remove(root + '.dem') 4149 os.remove(root + '.asc') 4150 os.remove(root + '.prj') 4151 4152 def test_dem2pts_bounding_box_Nullvalues(self): 4153 """Test conversion from dem in ascii format to native NetCDF xya format 4154 """ 4155 4156 import time, os 4157 from Numeric import array, zeros, allclose, Float, concatenate 4158 from Scientific.IO.NetCDF import NetCDFFile 4159 4160 #Write test asc file 4161 root = 'demtest' 4162 4163 filename = root+'.asc' 4164 fid = open(filename, 'w') 4165 fid.write("""ncols 5 4166 nrows 6 4167 xllcorner 2000.5 4168 yllcorner 3000.5 4169 cellsize 25 4170 NODATA_value 9999 4171 """) 4172 #Create linear function 4173 4174 ref_points = [] 4175 ref_elevation = [] 4176 new_ref_pts1 = [] 4177 new_ref_elev1 = [] 4178 NODATA_value = 9999 4179 for i in range(6): 4180 y = (6i)*25.0 4181 for j in range(5): 4182 x = j*25.0 4183 z = x+2*y 4184 if j == 4: z = NODATA_value # column 4185 if i == 2 and j == 2: z = NODATA_value # random 4186 if i == 5 and j == 1: z = NODATA_value 4187 if i == 1: z = NODATA_value # row 4188 if i == 3 and j == 1: z = NODATA_value # two pts/row 4189 if i == 3 and j == 3: z = NODATA_value 4190 4191 if z <> NODATA_value: 4192 new_ref_elev1.append(z) 4193 new_ref_pts1.append( [x,y] ) 4194 4195 ref_points.append( [x,y] ) 4196 ref_elevation.append(z) 4197 fid.write('%f ' %z) 4198 fid.write('\n') 4199 4200 fid.close() 4201 4202 #Write prj file with metadata 4203 metafilename = root+'.prj' 4204 fid = open(metafilename, 'w') 4205 4206 4207 fid.write("""Projection UTM 4208 Zone 56 4209 Datum WGS84 4210 Zunits NO 4211 Units METERS 4212 Spheroid WGS84 4213 Xshift 0.0000000000 4214 Yshift 10000000.0000000000 4215 Parameters 4216 """) 4217 fid.close() 4218 4219 #Convert to NetCDF pts 4220 convert_dem_from_ascii2netcdf(root) 4221 dem2pts(root, easting_min=2010.0, easting_max=2110.0, 4222 northing_min=3035.0, northing_max=3125.5) 4223 4224 #Check contents 4225 #Get NetCDF 4226 fid = NetCDFFile(root+'.pts', 'r') 4227 4228 # Get the variables 4229 #print fid.variables.keys() 4230 points = fid.variables['points'] 4231 elevation = fid.variables['elevation'] 4232 4233 #Check values 4234 assert fid.xllcorner[0] == 2010.0 4235 assert fid.yllcorner[0] == 3035.0 4236 4237 #create new reference points 4238 ref_points = [] 4239 ref_elevation = [] 4240 new_ref_pts2 = [] 4241 new_ref_elev2 = [] 4242 for i in range(4): 4243 y = (4i)*25.0 + 25.0 4244 y_new = y + 3000.5  3035.0 4245 for j in range(4): 4246 x = j*25.0 + 25.0 4247 x_new = x + 2000.5  2010.0 4248 z = x+2*y 4249 4250 if j == 3: z = NODATA_value # column 4251 if i == 1 and j == 1: z = NODATA_value # random 4252 if i == 4 and j == 0: z = NODATA_value 4253 if i == 0: z = NODATA_value # row 4254 if i == 2 and j == 0: z = NODATA_value # two pts/row 4255 if i == 2 and j == 2: z = NODATA_value 4256 4257 if z <> NODATA_value: 4258 new_ref_elev2.append(z) 4259 new_ref_pts2.append( [x_new,y_new] ) 4260 4261 4262 ref_points.append( [x_new,y_new] ) 4263 ref_elevation.append(z) 4264 4265 #print points[:] 4266 #print ref_points 4267 #assert allclose(points, ref_points) 4268 4269 #print attributes[:] 4270 #print ref_elevation 4271 #assert allclose(elevation, ref_elevation) 4272 4273 4274 assert len(points) == len(new_ref_pts2), 'length of returned points not correct' 4275 assert allclose(points, new_ref_pts2), 'points do not align' 4276 4277 #print 'elevation', elevation[:] 4278 assert len(elevation) == len(new_ref_elev2), 'length of returned elevation not correct' 4279 assert allclose(elevation, new_ref_elev2), 'elevations do not align' 4280 #Cleanup 4281 fid.close() 4282 4283 4284 os.remove(root + '.pts') 4285 os.remove(root + '.dem') 4286 os.remove(root + '.asc') 4287 os.remove(root + '.prj') 4288 4289 4290 4291 4292
Note: See TracChangeset
for help on using the changeset viewer.