[gdal-dev] Help with band pixel functions

[Stephen Woodbridge]

Hi Adriana,

Here are a few ideas:

import numpy as np

ima = gdal.Open(filea, 0)
banda = ima.GetRasterBand(1)
dataa = banda.ReadAsArray()

print type(dataa), dataa.shape, isinstance(dataa, np.ma.MaskedArray)

imb = gdal.Open(fileb, 0)
bandb = imb.GetRasterBand(1)
datab = bandb.ReadAsArray()

print type(datab), datab.shape, isinstance(datab, np.ma.MaskedArray)

So if you have masked arrays, where the fill_value=255, ie: your nodata 
value

The mask is True for nodata so

weights = dataa.mask | datab.mask
np.average([dataa, datab], axis=0, weights=weights)

You might need to invert the weights in which case try weights=~weights 
instead.

You might need to do this in three operations using appropriate masks
1) copy dataa values that are not masked but are in datab
2) copy datab values that are not masked but are in dataa
3) to copy the average of dataa and datab where they are both not masked.

Another option might be to use numpy.vectorize() to apply a function to 
each cell in an array:
https://docs.scipy.org/doc/numpy/reference/generated/numpy.vectorize.html

-Steve

On 6/24/2019 10:38 PM, Adriana Parra wrote:
> Hello,
>
> I have been trying to mosaic two images that partially overlap using 
> band pixel functions in Python. My images contain unsigned integers 
> (uint8) and the no data value is 255. I considered using the nanmean 
> function so that I would get average values in overlapping areas of 
> the two images and maintain the original value of the non-overlapping 
> areas, but for some reason the NAN values are being interpreted as 
> zero, resulting in half of the original value in non-overlapping 
> areas.  I would really appreciate any feedback. I am using GDAL 3.1.0 
> and python 3.6. Bellow is the code I used and I am sending attached 
> the input files.
> gdalbuildvrt  mosaic1.vrt -input_file_list  input_files.txt -srcnodata 
> "255"
>
>  import numpy as np
> def average(in_ar, out_ar, xoff, yoff, xsize, ysize, 
> raster_xsize,raster_ysize, buf_radius, gt, **kwargs):
>     np.round_(np.clip(np.nanmean(in_ar, axis = 0, dtype = 'uint8'),0,255),
>               out = out_ar)
>
> gdal_translate --config GDAL_VRT_ENABLE_PYTHON YES mosaic1.vrt mosaic1.tif
>
> Mosaic1.jpeg 
> <https://drive.google.com/a/nevada.unr.edu/file/d/1YgugcVnU6d5BmPYkYcotqBi3fFQcPpY8/view?usp=drive_web>
> raw_files.jpeg 
> <https://drive.google.com/a/nevada.unr.edu/file/d/1SQoH58V3MGtP-oXYnPYsf66snFOsWagC/view?usp=drive_web>
> test_files.7z 
> <https://drive.google.com/a/nevada.unr.edu/file/d/14No5nblxhNKeOU8k86Gfan7utYIUSqRU/view?usp=drive_web>
>
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev


---
This email has been checked for viruses by Avast antivirus software.
https://www.avast.com/antivirus