Re: Merge many raster using maximum value where overlaps

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

Re: Merge many raster using maximum value where overlaps

Rousseau Lambert2, Louis-Philippe (EC)
Hi everyone,

I would like to have your opinion on some processing I have to do.

So I have to merge many rasters (eventually more then 50). Some overlaps
(never completely), but others don't, my files when merged cover Canada.
I have to keep the maximum value where there is some overlap between
different files.

My first try was to use gdal_calc with --calc="maximum(A,B)", but it can
only input 26 files at a time. And The output of gdal_calc is only the
region of overlap, in my case I have to keep the whole extent of my
rasters, not only the region of overlap... Is there a way I could use
gdal_calc and keep everything, not only the overlap?

My second try was to use gdalwarp with a vrt of all my files as input
and use -r max for my resampling. My results were that it only
overlapped my geotiff and did not took the maximum value.

My last try, was to use gdalwarp to ensure that all my files have the
same extent and that pixel overlap perfectly. Then I create a mosaic
using gdal_merge of all my rasters. At this point I don't manage yet the
maximum value. Then I loop trough all my geotiff using gdal_calc to
output the overlap region between one geotiff and my mosaic. then I can
merge the overlap region with my mosaic. Finally I get a mosaic of the
maximum value where my raster overlap.

I don't really like this solution because if I have 50 files, I have to
use gdal_calc 50 times and merge the result 50 times... That seems a lot
of processing to simply get the maximum value in a mosaic.

Can you think of a better way to do? Any thoughts or opinions are welcome!

Louis-Philippe R. Lambert

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Merge many raster using maximum value where overlaps

Nicolas Cadieux
Hi,

You may be able to do that with the Saga grid mosaic.  You can select ´max’ for overlapping areas.  


 Nicolas

Le 17 nov. 2017 à 15:05, Rousseau Lambert2, Louis-Philippe (EC) <[hidden email]> a écrit :

Hi everyone,

I would like to have your opinion on some processing I have to do.

So I have to merge many rasters (eventually more then 50). Some overlaps
(never completely), but others don't, my files when merged cover Canada.
I have to keep the maximum value where there is some overlap between
different files.

My first try was to use gdal_calc with --calc="maximum(A,B)", but it can
only input 26 files at a time. And The output of gdal_calc is only the
region of overlap, in my case I have to keep the whole extent of my
rasters, not only the region of overlap... Is there a way I could use
gdal_calc and keep everything, not only the overlap?

My second try was to use gdalwarp with a vrt of all my files as input
and use -r max for my resampling. My results were that it only
overlapped my geotiff and did not took the maximum value.

My last try, was to use gdalwarp to ensure that all my files have the
same extent and that pixel overlap perfectly. Then I create a mosaic
using gdal_merge of all my rasters. At this point I don't manage yet the
maximum value. Then I loop trough all my geotiff using gdal_calc to
output the overlap region between one geotiff and my mosaic. then I can
merge the overlap region with my mosaic. Finally I get a mosaic of the
maximum value where my raster overlap.

I don't really like this solution because if I have 50 files, I have to
use gdal_calc 50 times and merge the result 50 times... That seems a lot
of processing to simply get the maximum value in a mosaic.

Can you think of a better way to do? Any thoughts or opinions are welcome!

Louis-Philippe R. Lambert

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Merge many raster using maximum value where overlaps

David Shean-2
Hi all,
This is also possible with the dem_mosaic tool, part of the NASA Ames Stereo Pipeline toolkit.  The utility is multithreaded and meant to handle an arbitrary number of large rasters efficiently.  You want the ‘—max’ command line option.

Unfortunately, the official NASA ASP site is currently down, but you can download precompiled binaries for Linux and OS X here: https://github.com/NeoGeographyToolkit/StereoPipeline/releases.  See the asp_book.pdf for documentation and usage on dem_mosaic.

Hope that helps.
-David


On Nov 20, 2017, at 1:48 PM, Nicolas Cadieux <[hidden email]> wrote:

Hi,

You may be able to do that with the Saga grid mosaic.  You can select ´max’ for overlapping areas.  


 Nicolas

Le 17 nov. 2017 à 15:05, Rousseau Lambert2, Louis-Philippe (EC) <[hidden email]> a écrit :

Hi everyone,

I would like to have your opinion on some processing I have to do.

So I have to merge many rasters (eventually more then 50). Some overlaps
(never completely), but others don't, my files when merged cover Canada.
I have to keep the maximum value where there is some overlap between
different files.

My first try was to use gdal_calc with --calc="maximum(A,B)", but it can
only input 26 files at a time. And The output of gdal_calc is only the
region of overlap, in my case I have to keep the whole extent of my
rasters, not only the region of overlap... Is there a way I could use
gdal_calc and keep everything, not only the overlap?

My second try was to use gdalwarp with a vrt of all my files as input
and use -r max for my resampling. My results were that it only
overlapped my geotiff and did not took the maximum value.

My last try, was to use gdalwarp to ensure that all my files have the
same extent and that pixel overlap perfectly. Then I create a mosaic
using gdal_merge of all my rasters. At this point I don't manage yet the
maximum value. Then I loop trough all my geotiff using gdal_calc to
output the overlap region between one geotiff and my mosaic. then I can
merge the overlap region with my mosaic. Finally I get a mosaic of the
maximum value where my raster overlap.

I don't really like this solution because if I have 50 files, I have to
use gdal_calc 50 times and merge the result 50 times... That seems a lot
of processing to simply get the maximum value in a mosaic.

Can you think of a better way to do? Any thoughts or opinions are welcome!

Louis-Philippe R. Lambert

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Merge many raster using maximum value where overlaps

Rutger
You could also consider using the VRT pixel functions (since 2.2). I'm not
sure about the performance when using the Python pixel functions (or perhaps
with Numba), i haven't done any benchmarks. It might be a 'nice to have' if
some statistics like max,min,std are added to the default functions if that
will be a lot faster. I'm not sure how much effort that would take, but
given that mul, diff etc are already there perhaps its really easy.

You can make a normal VRT mosaic with gdalbuildvrt, and then add the pixel
function to it yourself. The part you should add looks something like:



More information can be found at:
http://www.gdal.org/gdal_vrttut.html#gdal_vrttut_derived_python


Regards,
Rutger



--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Merge many raster using maximum value where overlaps

Rutger
Sorry, i had the code embedded as raw text. It appears to be left out of my
previous post, even though it worked in the preview.

Here it is:

  <VRTRasterBand dataType="Float32" band="1"
subClass="VRTDerivedRasterBand">
 
  <PixelFunctionType>max</PixelFunctionType>
  <PixelFunctionLanguage>Python</PixelFunctionLanguage>
  <PixelFunctionCode>
    </PixelFunctionCode>


Regards,
Rutger




--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Merge many raster using maximum value where overlaps

Rousseau Lambert2, Louis-Philippe (EC)
Thanks for the possible solution.

Unfortunately, I could not test using SAGA our NASA Ames Stereo Pipeline
toolkit. My script will have to be run on a operational environment in
which it is very difficult to have new tools installed.

Same thing for VRT pixel function. We don't have gdal 2.2 installed yet
(but it should be pretty soon) and will test this solution once I have
access to this version and let you know if it works for us!

For the moment the solution I have (different from the previous one and
way faster):

1- I reproject all my individual Geotiff with the same extent and
resolution and add a virtual value of 0 where I have no data.
2- Use a few iteration of gdal_calc (since we are limited to 26 files at
a time) to get the maximum value (now it ouputs all the extent of the
file because I populated it with 0 value and they all have the same extent).
3- Use a final gdal_calc to merge with the maximum value all the
previous gdal_calc reuslts.
4- Use gdal_translate to set 0 value pixels to no data (in my case, 0 is
really no data).

Thanks for your help!

Louis-Philippe R. Lambert

On 11/21/2017 09:34 AM, Rutger wrote:

> Sorry, i had the code embedded as raw text. It appears to be left out of my
> previous post, even though it worked in the preview.
>
> Here it is:
>
>   <VRTRasterBand dataType="Float32" band="1"
> subClass="VRTDerivedRasterBand">
>  
>   <PixelFunctionType>max</PixelFunctionType>
>   <PixelFunctionLanguage>Python</PixelFunctionLanguage>
>   <PixelFunctionCode>
>     </PixelFunctionCode>
>
>
> Regards,
> Rutger
>
>
>
>
> --
> Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
> _______________________________________________
> gdal-dev mailing list
> [hidden email]
> https://lists.osgeo.org/mailman/listinfo/gdal-dev


--
Louis-Philippe Rousseau Lambert, B.Sc.
Géomaticien / Geomatician
Section des Données, Performances et Standards
Data, Performance and Standards Section
Service Météorologique du Canada
Meteorological Service of Canada
Environnement et Changement Climatique Canada
Environment and Climate Change Canada
[hidden email]
(514) 421-5045

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev