Changeset 9553
- Timestamp:
- Jan 30, 2015, 9:22:42 AM (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/utilities/spatialInputUtil.py
r9513 r9553 725 725 rasterBand=raster.GetRasterBand(band) 726 726 rasterBandType=gdal.GetDataTypeName(rasterBand.DataType) 727 nodataval = rasterBand.GetNoDataValue() 727 728 728 729 # Projection info … … 788 789 else: 789 790 # Swap xl for xu 790 xu = xl 791 xu = xl + 0 791 792 xl = max(xu - 1, 0) 792 793 … … 794 795 yu = min(yl + 1, yMax - 1) 795 796 else: 796 yu = yl 797 yu = yl + 0 797 798 yl = max(yu - 1, 0) 798 799 799 800 # Map x,y to unit square 800 x = px[i] - (xl + 0.5) 801 y = py[i] - (yl + 0.5) 801 if(xu > xl): 802 x = px[i] - (xl + 0.5) 803 else: 804 x = 0. 805 806 if(yu > yl): 807 y = py[i] - (yl + 0.5) 808 else: 809 y = 0. 802 810 803 811 if not ( (x>=0.) & (x<=1.)): 804 print x, xl, xu, px[i]812 print 'x-values error: ', x, xl, xu, px[i], xMax 805 813 raise Exception('x out of bounds') 806 814 807 815 if not ( (y>=0.) & (y<=1.)): 808 print y, yl, yu, py[i]816 print 'y-values error: ', y, yl, yu, py[i] 809 817 raise Exception('y out of bounds') 810 818 … … 827 835 828 836 # Bilinear interpolation 829 elev[i] = r00*(1-x)*(1-y) + r01*(1-x)*y +\ 830 r10*x*(1-y) + r11*x*y 837 elev[i] = r00*(1.-x)*(1.-y) + r01*(1.-x)*y +\ 838 r10*x*(1.-y) + r11*x*y 839 840 # Deal with nodata. This needs to be in the loop 841 # Just check if any of the pixels are nodata 842 if nodataval is not None: 843 if numpy.isfinite(nodataval): 844 rij = numpy.array([r00, r01, r10, r11]) 845 rel_tol = ( abs(rij - nodataval) < nodata_rel_tol*abs(nodataval) ) 846 missing = (rel_tol).nonzero()[0] 847 if len(missing) > 0: 848 elev[i] = numpy.nan 831 849 else: 832 850 raise Exception('Unknown value of "interpolation"') 833 851 834 # Deal with nodata 835 nodataval = rasterBand.GetNoDataValue() 836 if nodataval is not None: 837 if numpy.isfinite(nodataval): 838 rel_tol = ( abs(elev - nodataval) < nodata_rel_tol*abs(nodataval) ) 839 missing = (rel_tol).nonzero()[0] 840 if len(missing) > 0: 841 elev[missing] = numpy.nan 852 # Deal with nodata for pixel based interpolation [efficient treatment 853 # outside of loop] 854 if (interpolation == 'pixel'): 855 if nodataval is not None: 856 if numpy.isfinite(nodataval): 857 rel_tol = ( abs(elev - nodataval) < nodata_rel_tol*abs(nodataval) ) 858 missing = (rel_tol).nonzero()[0] 859 if len(missing) > 0: 860 elev[missing] = numpy.nan 842 861 843 862 return elev
Note: See TracChangeset
for help on using the changeset viewer.