Changeset 9553


Ignore:
Timestamp:
Jan 30, 2015, 9:22:42 AM (9 years ago)
Author:
davies
Message:

Bugfix related to bilinear raster lookup

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/utilities/spatialInputUtil.py

    r9513 r9553  
    725725        rasterBand=raster.GetRasterBand(band)
    726726        rasterBandType=gdal.GetDataTypeName(rasterBand.DataType)
     727        nodataval = rasterBand.GetNoDataValue()
    727728   
    728729        # Projection info
     
    788789                else:
    789790                    # Swap xl for xu
    790                     xu = xl
     791                    xu = xl + 0
    791792                    xl = max(xu - 1, 0)
    792793
     
    794795                    yu = min(yl + 1, yMax - 1)
    795796                else:
    796                     yu = yl
     797                    yu = yl + 0
    797798                    yl = max(yu - 1, 0)
    798799
    799800                # 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.
    802810
    803811                if not ( (x>=0.) & (x<=1.)):
    804                     print x, xl, xu, px[i]
     812                    print 'x-values error: ', x, xl, xu, px[i], xMax
    805813                    raise Exception('x out of bounds')
    806814
    807815                if not ( (y>=0.) & (y<=1.)):
    808                     print y, yl, yu, py[i]
     816                    print 'y-values error: ', y, yl, yu, py[i]
    809817                    raise Exception('y out of bounds')
    810818
     
    827835
    828836                # 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
    831849            else:
    832850                raise Exception('Unknown value of "interpolation"')           
    833851
    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
    842861
    843862        return elev
Note: See TracChangeset for help on using the changeset viewer.