search

OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.91k stars 2.55k forks source link

how to blend two dem? #11276

Open sendreams opened 6 hours ago

sendreams commented 6 hours ago

Feature description

i have two dem tif files which should be blend to a new tif file,

the files information is: file1 (no compress)

C:\Apps\dev\gdal\release-1928-x64-gdal-3-9-1-mapserver-8-0-1\bin\gdal\apps>gdalinfo C:\data\dem\blend\1030-UTM49.tif
Driver: GTiff/GeoTIFF
Files: C:\data\dem\blend\1030-UTM49.tif
       C:\data\dem\blend\1030-UTM49.tif.ovr
       C:\data\dem\blend\1030-UTM49.tif.aux.xml
Size is 46657, 24600
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (109.163925414663652,29.539997383024080)
Pixel Size = (0.000004511081671,-0.000004511081671)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 109.1639254,  29.5399974) (109d 9'50.13"E, 29d32'23.99"N)
Lower Left  ( 109.1639254,  29.4290248) (109d 9'50.13"E, 29d25'44.49"N)
Upper Right ( 109.3743990,  29.5399974) (109d22'27.84"E, 29d32'23.99"N)
Lower Right ( 109.3743990,  29.4290248) (109d22'27.84"E, 29d25'44.49"N)
Center      ( 109.2691622,  29.4845111) (109d16' 8.98"E, 29d29' 4.24"N)
Band 1 Block=46657x1 Type=Float32, ColorInterp=Gray
  Min=414.073 Max=952.708
  Minimum=414.073, Maximum=952.708, Mean=667.534, StdDev=147.463
  NoData Value=-32767
  Overviews: 23329x12300, 11665x6150, 5833x3075, 2917x1538, 1459x769, 730x385, 365x193, 183x97, 92x49
  Unit Type: metre
  Metadata:
    STATISTICS_MAXIMUM=952.70776367188
    STATISTICS_MEAN=667.53419949254
    STATISTICS_MINIMUM=414.07342529297
    STATISTICS_STDDEV=147.46325087621
    STATISTICS_VALID_PERCENT=1.309

file2 (COMPRESSION=LZW)

C:\Apps\dev\gdal\release-1928-x64-gdal-3-9-1-mapserver-8-0-1\bin\gdal\apps>gdalinfo C:\data\dem\blend\TestRegion.tif
Driver: GTiff/GeoTIFF
Files: C:\data\dem\blend\TestRegion.tif
       C:\data\dem\blend\TestRegion.tif.ovr
       C:\data\dem\blend\TestRegion.tif.xml
       C:\data\dem\blend\TestRegion.tif.aux.xml
Size is 5072, 2979
Coordinate System is:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (109.161255799295816,29.549341546442445)
Pixel Size = (0.000045105627232,-0.000045105964498)
Metadata:
  AREA_OR_POINT=Area
  DataType=Generic
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 109.1612558,  29.5493415) (109d 9'40.52"E, 29d32'57.63"N)
Lower Left  ( 109.1612558,  29.4149709) (109d 9'40.52"E, 29d24'53.90"N)
Upper Right ( 109.3900315,  29.5493415) (109d23'24.11"E, 29d32'57.63"N)
Lower Right ( 109.3900315,  29.4149709) (109d23'24.11"E, 29d24'53.90"N)
Center      ( 109.2756437,  29.4821562) (109d16'32.32"E, 29d28'55.76"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
  Min=409.893 Max=1019.884
  Minimum=409.893, Maximum=1019.884, Mean=619.570, StdDev=115.152
  NoData Value=-3.402823e+38
  Overviews: 2536x1490, 1268x745, 634x373, 317x187
  Metadata:
    RepresentationType=ATHEMATIC
    STATISTICS_COVARIANCES=13259.90891525758
    STATISTICS_MAXIMUM=1019.8840332031
    STATISTICS_MEAN=619.5701603997
    STATISTICS_MEDIAN=598.870576
    STATISTICS_MINIMUM=409.89294433594
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=115.15167786558

blend two files, use "COMPRESS=LZW" parameter

gdalwarp -t_srs epsg:4326 -r cubicspline -co COMPRESS=LZW -co BIGTIFF=YES -co TILED=YES -co BLOCKXSIZE=128 -co BLOCKYSIZE=128 C:/data/dem/blend/TestRegion.tif C:/data/dem/blend/1030-UTM49.tif C:/data/dem/blend/blend1.tif

the output tif info: ( COMPRESSION=LZW)

C:\Apps\dev\gdal\release-1928-x64-gdal-3-9-1-mapserver-8-0-1\bin\gdal\apps>gdalinfo C:\data\dem\blend\blend1.tif
Driver: GTiff/GeoTIFF
Files: C:\data\dem\blend\blend1.tif
Size is 50714, 29787
Coordinate System is:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (109.161255799295816,29.549341546442445)
Pixel Size = (0.000004511081671,-0.000004511081671)
Metadata:
  DataType=Generic
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 109.1612558,  29.5493415) (109d 9'40.52"E, 29d32'57.63"N)
Lower Left  ( 109.1612558,  29.4149700) (109d 9'40.52"E, 29d24'53.89"N)
Upper Right ( 109.3900308,  29.5493415) (109d23'24.11"E, 29d32'57.63"N)
Lower Right ( 109.3900308,  29.4149700) (109d23'24.11"E, 29d24'53.89"N)
Center      ( 109.2756433,  29.4821558) (109d16'32.32"E, 29d28'55.76"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
  NoData Value=-3.402823e+38
  Metadata:
    RepresentationType=ATHEMATIC

all files size:

1731601801918

then i try to using ArcGIS Pro to blend the same tif files

1731602068441

the result tif info : ( COMPRESSION=LZW)

C:\Apps\dev\gdal\release-1928-x64-gdal-3-9-1-mapserver-8-0-1\bin\gdal\apps>gdalinfo C:\data\dem\blend\mosaic_arcgis.tif
Driver: GTiff/GeoTIFF
Files: C:\data\dem\blend\mosaic_arcgis.tif
       C:\data\dem\blend\mosaic_arcgis.tif.ovr
       C:\data\dem\blend\mosaic_arcgis.tif.xml
       C:\data\dem\blend\mosaic_arcgis.tif.aux.xml
Size is 50715, 29788
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (109.161254854314251,29.549344344246968)
Pixel Size = (0.000004511081671,-0.000004511081671)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 109.1612549,  29.5493443) (109d 9'40.52"E, 29d32'57.64"N)
Lower Left  ( 109.1612549,  29.4149682) (109d 9'40.52"E, 29d24'53.89"N)
Upper Right ( 109.3900344,  29.5493443) (109d23'24.12"E, 29d32'57.64"N)
Lower Right ( 109.3900344,  29.4149682) (109d23'24.12"E, 29d24'53.89"N)
Center      ( 109.2756446,  29.4821563) (109d16'32.32"E, 29d28'55.76"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
  Min=409.893 Max=1019.884
  Minimum=409.893, Maximum=1019.884, Mean=619.440, StdDev=115.091
  NoData Value=3.4e+38
  Overviews: 25358x14894, 12679x7447, 6340x3724, 3170x1862, 1585x931, 793x466, 397x233
  Unit Type: metre
  Metadata:
    RepresentationType=ATHEMATIC
    STATISTICS_COVARIANCES=13246.0194951894
    STATISTICS_MAXIMUM=1019.8840332031
    STATISTICS_MEAN=619.44029421632
    STATISTICS_MEDIAN=596.478454
    STATISTICS_MINIMUM=409.89294433594
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=115.09135282544

these two files information almost identical, but the arcgis output file size is much smaller.

1731602425751

what's the reason about it?

Additional context

gdal version is : 3.9.1

rouault commented 6 hours ago

gdalwarp processes files one by one, so if they overlap significantly, it may end up updating the same output tile several times, which in case of a compressed dataset, leads to space being lost in the file

You should probably use: gdalbuildvrt -r cubic tmp.vrt input1.tif input2.tif followed by gdal_translate tmp.vrt out.tif -co TILED=YES -co COMPRESS=LZW

sendreams commented 5 hours ago

@rouault thank you so much, you save my life.

rouault commented 5 hours ago

as the workflow of gdalbuildvrt + gdal_translate is not totally equilvalent as the one of gdalwarp, if you need to use gdalwarp, one alternative is just to add a post-processing stage with gdal_translate output_of_gdalwarp.tif final_result.tif -co TILED=YES -co COMPRESS=LZW that will generate a compact raster with the same content as the one output by gdalwarp