Positive and negative masks

I'm working on kriging at the moment and need to split my beloved glacier into two regions: the interpolated region and the extrapolated region.

1. Firstly, create a delaunay polygon shape file in QGIS from point data. This represents the area of the glacier where values can be interpolated (with some editing). 
2. From this and a mask of the entire glacier (a raster with 1 = glacier, nan = not a glacier) I can use gdalwarp to create a masked mask, i.e. mask the glacier to show only those parts covered by interpolation.
3.  gdalwarp refuses to respect the size of the original raster when clipping, so that has to be fixed running gdalwarp again and forcing the new raster to the right size
4. Now, to get the reverse mask, the bits not covered by interpolation. Well, there I hit a bit of a problem. gdal_calc.py won't do it, at least no matter what I tried.
5. So write a script.

Maybe I'm missing something but by this point, writing something myself seemed quicker and infinitely less frustrating.

The function below uses a couple of other functions already encountered here on this blog, the DemImport function and namer, although you can easily get rid of namer. The code is messy and probably not efficient but I'm bolting together this stuff on the fly, as needs arise. But then, it's free.

def maskPosNeg(vector,Afile,name,outfolder):
    '''Create a mask raster from a vector over a raster then invert masked area to create mask of remaining area
    Created for masks of interpolated and extrapolated ares of glacier
    To use: maskPosNeg(vector,Afile,name,outfolder) '''
    # Get background raster (glacier dem) and calculate dimensions etc.
    Adata,AXg,AYg,AXg1,AYg1,Arx,Ary,Ademmeta = DemImport(Afile)
    transf = Ademmeta[2]
    ul_x = transf[0]
    ul_y = transf[3]
    xres = transf[1]
    yres = transf[5]
    # get image size
    demrows = Ademmeta[4]
    demcols = Ademmeta[5]
    # Calculate corners
    ll_x = ul_x
    ll_y = ul_y + (demrows * yres)
    ur_x = ll_x + (demcols * xres)
    ur_y = ul_y
    # Create names for output data
    Aname,Aext,Apath,Anamefull = namer(Afile)
    outFile1 = os.path.join(outfolder,(name + '_posMask_Temp.tif'))
    outFile2 = os.path.join(outfolder,(name + '_posMask.tif'))
    outFile3 = name + '_negMask'
    # Start by clipping to vector
    gdalMess1 = 'gdalwarp -te ' + str(ll_x) + ' ' + str(ll_y) + ' ' + str(ur_x) + ' ' + str(ur_y) + ' -tr ' +str(xres) + ' ' + str(yres) + ' -dstnodata -9999 -q -cutline ' + vector + ' -crop_to_cutline -of GTiff ' + Afile + ' ' + outFile1
    os.system(gdalMess1)
    # Then resize to original raster, I kid you not.
    gdalMess2 = 'gdalwarp -te ' + str(ll_x) + ' ' + str(ll_y) + ' ' + str(ur_x) + ' ' + str(ur_y) + ' -tr ' +str(xres) + ' ' + str(yres) + ' -dstnodata -9999 ' + outFile1 + ' ' + outFile2
    os.system(gdalMess2)
    print outFile2,' created'
    # Remove temp file
    command = 'rm ' + outFile1
    os.system(command)
    #
    Bfile = outFile2
    Bdata,BXg,BYg,BXg1,BYg1,Brx,Bry,Bdemmeta = DemImport(Bfile)
    print Bfile,' data type: ',type(Bdata)
    print 'Rows (A,B): ',Arx,Brx,' Columns (A,B): ',Ary,Bry
    print 'xres (A,B): ',Ademmeta[2][1],Bdemmeta[2][1],' yres (A,B): ',Ademmeta[2][5],Bdemmeta[2][5]
    print 'No data (A): ',Ademmeta[6],' No data (B): ',Bdemmeta[6]
    # Check matching resolution
    if Arx != Brx or Ary != Bry:
        print Bfile,' resolution mismatch with ', Afile
    elif Ademmeta[4] != Bdemmeta[4] or Ademmeta[5] != Bdemmeta[5]:
        print 'Size mismatch between ',Bfile,' and ',Afile
    # Do the masking and subtracting
    Bdata[Bdata==Bdemmeta[6]]=np.nan
    AdataMasked = np.ma.masked_array(np.nan_to_num(Adata), mask=np.isnan(Adata))# & np.isnan(Bdata))
    BdataMasked = np.ma.masked_array(np.nan_to_num(Bdata), mask=AdataMasked.mask)
    outdata = (AdataMasked - BdataMasked).filled(np.nan)
    outdata = np.ma.masked_array(np.nan_to_num(outdata), mask=np.isnan(outdata))
    filnm = datawrite(outdata,Adata,Ademmeta,outFile3,outfolder)
    print 'maskPosNeg done.'
    return 0