[gdal-dev] gdalbuildvrt Max Value / Pixel Function

2025-02-13 Thread Abdul Raheem Siddiqui via gdal-dev
Hello,

I am looking for functionality to create a VRT of multiple rasters where
the VRT pixel value is the maximum of underlying rasters pixel values.

I found out that we have pixel functions
https://lists.osgeo.org/pipermail/gdal-dev/2016-September/045134.html that
we can add to VRTs to get this behavior but just checking to see if there
is a direct way to specify this in gdalbuildvrt command line call or do one
has to manually inject this in VRT xml.

On the same note does the new gdaltileindex provide this functionality of
selecting the max value for each pixel out of the box?

Regards,
Abdul Raheem Siddiqui
ERT Inc.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev


Re: [gdal-dev] Debugging mask copy issues when using Band interleaved COGs

2025-02-13 Thread Even Rouault via gdal-dev

Owen,

Please create a ticket about that.

Even
Le 13/02/2025 à 04:40, Owen Glofcheski via gdal-dev a écrit :

Hey all,

I work with Kurt Schwehr, and we've been debugging failures where an 
explicit mask isn't preserved through a CreateCopy operation to create 
a COG through the GTiff driver at HEAD.


This failure only started occurring after 
https://github.com/OSGeo/gdal/commit/3bdc81a205bad8342688873e4ed5716f1208e8f5.


After the CreateCopy, when we read back the mask, the mask is expanded 
from the 1-bit mask to the 8-bit mask inside of gtiffoddbitsband.h. 
However, the mask is not a 1-bit mask when I add logging here; 
instead, it's the original 8-bit mask (with 255s and 0s). So when it 
gets expanded to 8-bits, it's no longer the original mask:

E.g., 255 0 255 255 becomes 255 255 ... 255 255 0 0 ... 0 0 255 255 ...

So my guess is that somewhere along the line we either copy the mask 
without turning it back into a bitmask, or we lose metadata tracking 
that the mask is encoded a certain way.


Does anything obvious jump out to folks here? The issue *does not* 
occur when either COPY_SRC_OVERVIEWS="YES" is removed or INTERLEAVE is 
set to PIXEL (it only fails on band interleaving).


Here's a gist containing an example test that fails at head but 
previously succeeded: 
https://gist.github.com/orglofch/40cf4d9136978b0384141cc74d1039b6.


Thanks,
Owen

___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev


--
http://www.spatialys.com
My software is free, but my time generally not.

___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev


Re: [gdal-dev] Problem with warping into existing TIF

2025-02-13 Thread Even Rouault via gdal-dev

Tim,

Reviewing the warping logic, I see we have support for a *input* mask, 
but not creating/updating an *output* mask. AFAICS, there isn't any 
particular strong reason not to have support for that. Likely just a few 
hours of coding away. Your best workaround is to generate an output with 
an alpha band (possibly as a VRT) and use gdal.Translate(dst, src, 
options="-b 1 -b 2 -b 3 -mask 4") for example to transform a RGBA into a 
RGB+mask product


Even

Le 12/02/2025 à 23:47, Tim Harris via gdal-dev a écrit :
Not sure if this is a bug, or expected behavior, or user error. If I 
warp one TIF into another, and both have nodata masks, it seems that 
gdalwarp isn't updating the destination TIF's nodata mask to unmask 
the new pixels.


Here's a short script to generate two example TIFs. One has a red box 
in its upper left corner, the other has a green box in its upper right 
corner. Both have nodata masks so that only the colored boxes are visible.


---
import numpy as np
from osgeo import gdal

gdal.UseExceptions()
gdal.SetConfigOption("GDAL_TIFF_INTERNAL_MASK", "YES")
drv = gdal.GetDriverByName("GTiff")

zeros = np.zeros((3, 1024, 1024))

red = np.copy(zeros)
red[0, 100:200, 100:200] = 255

green = np.copy(zeros)
green[1, 100:200, 800:900] = 255

ds = drv.Create("red.tif", 1024, 1024, 3, gdal.GDT_Byte)
ds.SetGeoTransform((0, 1, 0, 0, 0, -1))
ds.WriteArray(red)
ds.CreateMaskBand(gdal.GMF_PER_DATASET)
ds.GetRasterBand(1).GetMaskBand().WriteArray(red[0, :, :])
ds = None

ds = drv.Create("green.tif", 1024, 1024, 3, gdal.GDT_Byte)
ds.SetGeoTransform((0, 1, 0, 0, 0, -1))
ds.WriteArray(green)
ds.CreateMaskBand(gdal.GMF_PER_DATASET)
ds.GetRasterBand(1).GetMaskBand().WriteArray(green[1, :, :])
ds = None
---

Then if you copy red.tif to warp.tif, then warp green.tif into warp.tif:
cp red.tif warp.tif
gdalwarp green.tif warp.tif

I would expect the result to be a combined image with both the red and 
green boxes. But when I load this result into QGIS, I only see the red 
box. If I change the layer in QGIS to ignore the mask band, I do see 
both boxes (and the remaining image is black, as expected). So the 
pixel data is there, but the mask didn't get updated.


Is there a way to get gdalwarp to update the nodata mask 
automatically? Maybe some hidden setting I'm missing? I can work 
around this by separately warping the masks, then copying the result 
into the final TIF's mask band. But it would be nice if the warp took 
care of this for me.


Thanks

___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev


--
http://www.spatialys.com
My software is free, but my time generally not.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev


Re: [gdal-dev] gdalbuildvrt Max Value / Pixel Function

2025-02-13 Thread Daniel Baston via gdal-dev
>
> Any preliminary documentation on this?
>

Preliminary, yes:
https://gdal--11850.org.readthedocs.build/en/11850/programs/gdal_raster_calc.html
(from https://github.com/OSGeo/gdal/pull/11850)

Dan
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev


Re: [gdal-dev] gdalbuildvrt Max Value / Pixel Function

2025-02-13 Thread Scott via gdal-dev

WOW! Game changer.

Any preliminary documentation on this?


On 2/13/25 13:00, Daniel Baston via gdal-dev wrote:
the next release of GDAL should have a replacement "gdal_calc" that 
would take multiple inputs and write a VRT that applies an expression to 
those inputs


___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev


Re: [gdal-dev] gdalbuildvrt Max Value / Pixel Function

2025-02-13 Thread Daniel Baston via gdal-dev
Hi Abdul,

For the specific case of multiple rasters that share the same extent, the
next release of GDAL should have a replacement "gdal_calc" that would take
multiple inputs and write a VRT that applies an expression to those inputs
("max", in your case). In the meantime, I think your best option would be
to prepare the VRT yourself.

Dan

On Thu, Feb 13, 2025 at 3:12 PM Abdul Raheem Siddiqui via gdal-dev <
gdal-dev@lists.osgeo.org> wrote:

> Hello,
>
> I am looking for functionality to create a VRT of multiple rasters where
> the VRT pixel value is the maximum of underlying rasters pixel values.
>
> I found out that we have pixel functions
> https://lists.osgeo.org/pipermail/gdal-dev/2016-September/045134.html
> that we can add to VRTs to get this behavior but just checking to see if
> there is a direct way to specify this in gdalbuildvrt command line call or
> do one has to manually inject this in VRT xml.
>
> On the same note does the new gdaltileindex provide this functionality of
> selecting the max value for each pixel out of the box?
>
> Regards,
> Abdul Raheem Siddiqui
> ERT Inc.
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev


Re: [gdal-dev] gdalbuildvrt Max Value / Pixel Function

2025-02-13 Thread Abdul Raheem Siddiqui via gdal-dev
Thanks for your reply, my rasters are not sharing extent but some do
overlap.

How about GTI, Does that provide this functionality? Do you know.

On Thu, Feb 13, 2025, 4:00 PM Daniel Baston  wrote:

> Hi Abdul,
>
> For the specific case of multiple rasters that share the same extent, the
> next release of GDAL should have a replacement "gdal_calc" that would take
> multiple inputs and write a VRT that applies an expression to those inputs
> ("max", in your case). In the meantime, I think your best option would be
> to prepare the VRT yourself.
>
> Dan
>
> On Thu, Feb 13, 2025 at 3:12 PM Abdul Raheem Siddiqui via gdal-dev <
> gdal-dev@lists.osgeo.org> wrote:
>
>> Hello,
>>
>> I am looking for functionality to create a VRT of multiple rasters where
>> the VRT pixel value is the maximum of underlying rasters pixel values.
>>
>> I found out that we have pixel functions
>> https://lists.osgeo.org/pipermail/gdal-dev/2016-September/045134.html
>> that we can add to VRTs to get this behavior but just checking to see if
>> there is a direct way to specify this in gdalbuildvrt command line call or
>> do one has to manually inject this in VRT xml.
>>
>> On the same note does the new gdaltileindex provide this functionality of
>> selecting the max value for each pixel out of the box?
>>
>> Regards,
>> Abdul Raheem Siddiqui
>> ERT Inc.
>>
>> ___
>> gdal-dev mailing list
>> gdal-dev@lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev