Changeset 2555
 Timestamp:
 Mar 17, 2006, 2:03:19 PM (18 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

inundation/pyvolution/test_data_manager.py
r2553 r2555 651 651 #create new reference points 652 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)653 newz[0:5] = ref_elevation[32:38] 654 newz[6:11] = ref_elevation[42:48] 655 newz[12:17] = ref_elevation[52:58] 656 newz[18:23] = ref_elevation[62:68] 657 657 ref_elevation = [] 658 658 ref_elevation = newz 659 659 ref_points = [] 660 new_ref_points = []661 660 x0 = 2002 662 661 y = 3007 … … 669 668 x = x0 + xvec[j] 670 669 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) 679 680 #print attributes[:] 681 #print 'ref_elev', ref_elevation[:] 682 #print 'received elev', elevation[:] 670 ref_points.append ([xnew,ynew]) #Relative point values 671 672 assert allclose(points, ref_points) 673 683 674 assert allclose(elevation, ref_elevation) 684 675 … … 720 711 yvec = range(10) 721 712 xvec = range(10) 722 z = range(100) 713 #z = range(100) 714 z = zeros(100) 723 715 NODATA_value = 9999 724 716 count = 1 … … 729 721 ref_points.append ([x,y]) 730 722 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 723 z[count] = (4*i  3*j)%13 724 if j == 4: z[count] = NODATA_value #column inside clipping region 725 if j == 8: z[count] = NODATA_value #column outside clipping region 726 if i == 9: z[count] = NODATA_value #row outside clipping region 727 if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region 734 728 ref_elevation.append( z[count] ) 735 729 fid.write('%f ' %z[count]) … … 776 770 777 771 #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 = [] 772 newz = zeros(19) 773 newz[0:2] = ref_elevation[32:34] 774 newz[2:5] = ref_elevation[35:38] 775 newz[5:7] = ref_elevation[42:44] 776 newz[7] = ref_elevation[45] 777 newz[8] = ref_elevation[47] 778 newz[9:11] = ref_elevation[52:54] 779 newz[11:14] = ref_elevation[55:58] 780 newz[14:16] = ref_elevation[62:64] 781 newz[16:19] = ref_elevation[65:68] 782 783 789 784 ref_elevation = newz 790 785 ref_points = [] … … 800 795 x = x0 + xvec[j] 801 796 xnew = x  2002.0 802 if j <> 2 :797 if j <> 2 and (i<>1 or j<>4): 803 798 ref_points.append([x,y]) 804 new_ref_points.append ([xnew,ynew]) 799 new_ref_points.append ([xnew,ynew]) 800 805 801 806 802 assert allclose(points, new_ref_points) 807 803 assert allclose(elevation, ref_elevation) 804 805 #Cleanup 806 fid.close() 807 808 809 os.remove(root + '.pts') 810 os.remove(root + '.dem') 811 os.remove(root + '.asc') 812 os.remove(root + '.prj') 813 814 815 def test_dem2pts_bounding_box_removeNullvalues_v3(self): 816 """Test conversion from dem in ascii format to native NetCDF xya format 817 Check missing values on clipping boundary 818 """ 819 820 import time, os 821 from Numeric import array, zeros, allclose, Float, concatenate, ones 822 from Scientific.IO.NetCDF import NetCDFFile 823 824 #Write test asc file 825 root = 'demtest' 826 827 filename = root+'.asc' 828 fid = open(filename, 'w') 829 fid.write("""ncols 10 830 nrows 10 831 xllcorner 2000 832 yllcorner 3000 833 cellsize 1 834 NODATA_value 9999 835 """) 836 #Create linear function 837 ref_points = [] 838 ref_elevation = [] 839 x0 = 2000 840 y = 3010 841 yvec = range(10) 842 xvec = range(10) 843 #z = range(100) 844 z = zeros(100) 845 NODATA_value = 9999 846 count = 1 847 for i in range(10): 848 y = y  1 849 for j in range(10): 850 x = x0 + xvec[j] 851 ref_points.append ([x,y]) 852 count += 1 853 z[count] = (4*i  3*j)%13 854 if j == 4: z[count] = NODATA_value #column inside clipping region 855 if j == 8: z[count] = NODATA_value #column outside clipping region 856 if i == 6: z[count] = NODATA_value #row on clipping boundary 857 if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region 858 ref_elevation.append( z[count] ) 859 fid.write('%f ' %z[count]) 860 fid.write('\n') 861 862 fid.close() 863 864 #print 'sending elev', ref_elevation 865 866 #Write prj file with metadata 867 metafilename = root+'.prj' 868 fid = open(metafilename, 'w') 869 870 871 fid.write("""Projection UTM 872 Zone 56 873 Datum WGS84 874 Zunits NO 875 Units METERS 876 Spheroid WGS84 877 Xshift 0.0000000000 878 Yshift 10000000.0000000000 879 Parameters 880 """) 881 fid.close() 882 883 #Convert to NetCDF pts 884 convert_dem_from_ascii2netcdf(root) 885 dem2pts(root, easting_min=2002.0, easting_max=2007.0, 886 northing_min=3003.0, northing_max=3006.0) 887 888 #Check contents 889 #Get NetCDF 890 fid = NetCDFFile(root+'.pts', 'r') 891 892 # Get the variables 893 #print fid.variables.keys() 894 points = fid.variables['points'] 895 elevation = fid.variables['elevation'] 896 897 #Check values 898 assert fid.xllcorner[0] == 2002.0 899 assert fid.yllcorner[0] == 3003.0 900 901 #create new reference points 902 newz = zeros(14) 903 newz[0:2] = ref_elevation[32:34] 904 newz[2:5] = ref_elevation[35:38] 905 newz[5:7] = ref_elevation[42:44] 906 newz[7] = ref_elevation[45] 907 newz[8] = ref_elevation[47] 908 newz[9:11] = ref_elevation[52:54] 909 newz[11:14] = ref_elevation[55:58] 910 911 912 913 ref_elevation = newz 914 ref_points = [] 915 new_ref_points = [] 916 x0 = 2002 917 y = 3007 918 yvec = range(4) 919 xvec = range(6) 920 for i in range(4): 921 y = y  1 922 ynew = y  3003.0 923 for j in range(6): 924 x = x0 + xvec[j] 925 xnew = x  2002.0 926 if j <> 2 and (i<>1 or j<>4) and i<>3: 927 ref_points.append([x,y]) 928 new_ref_points.append ([xnew,ynew]) 929 930 931 #print points[:],points[:].shape 932 #print new_ref_points, len(new_ref_points) 933 934 assert allclose(elevation, ref_elevation) 935 assert allclose(points, new_ref_points) 936 808 937 809 938 #Cleanup
Note: See TracChangeset
for help on using the changeset viewer.