3

I'm trying to run the following code to downsample a ~3GB GTiff while coregistering it to a reference extent:

def coreg_downsamp(source, destination, output_type="GTiff"):
    #essentially we're Warping into the same CRS, but with a resolution of 10m,
    # and with destination extent defined by the ELME reference raster
    # (also, using multi-threaded processing and finally compression)
    
    #attempt with python bindings
    #opts =  f"-tr 10.0 10.0 -r med -te {REFERENCE_EXTENT_STR} -dstnodata -9999.0 -multi -of {output_type} -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9"
    #gdal.Warp(destNameOrDestDS = destination, srcDSOrSrcDSTab = source, options=opts)
    
    #just sending it straight to command line
    command = f"gdalwarp -dstnodata -9999.0 -tr 10.0 10.0 -r med -te {REFERENCE_EXTENT_STR} -multi -of {output_type} -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 {source} {destination}"
    subprocess.call(command, shell=True)

Raw command: gdalwarp -dstnodata -9999.0 -tr 10.0 10.0 -r med -te 369000.0 6377000.0 740000.0 6635000.0 -multi -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 eesti_ndvi_suvi2022.tif ../processed_data/ndvitest_commandline.tif

I built the command using the QGIS Warp tool, which you can see below:

enter image description here

Running the QGIS GDAL Warp tool works great! It does everything I need, and the result is even about ~2GB after compression. However, if I take the exact same command and run it in the command line, send it to the command line from Python, or use the Python GDAL bindings, the result is more like ~10GB and I have to turn on the BIGTIFF flag for it to even work: enter image description here

The only difference at all is that I'm running my command line/Python in Linux (WSL2 Ubuntu) whereas my QGIS is running on Windows, so I've changed the I/O filepaths accordingly. Both are using GDAL 3.4.3.

For good measure, here's the GDALInfo on the input file:

Driver: GTiff/GeoTIFF
Files: eesti_ndvi_suvi2022.tif
Size is 37013, 24284
Coordinate System is:
PROJCRS["Estonian Coordinate System of 1997",
    BASEGEOGCRS["EST97",
        DATUM["Estonia 1997",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4180]],
    CONVERSION["Estonian National Grid",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",57.5175539305556,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",24,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",59.3333333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",58,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",6375000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["northing (X)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (Y)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Topographic mapping (large scale)."],
        AREA["Estonia - onshore and offshore."],
        BBOX[57.52,20.37,60,28.2]],
    ID["EPSG",3301]]
Data axis to CRS axis mapping: 2,1
Origin = (369010.000000000000000,6620000.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  369010.000, 6620000.000) ( 21d40'24.68"E, 59d41'48.08"N)
Lower Left  (  369010.000, 6377160.000) ( 21d48'50.50"E, 57d31' 4.02"N)
Upper Right (  739140.000, 6620000.000) ( 28d14'36.00"E, 59d38'57.00"N)
Lower Right (  739140.000, 6377160.000) ( 27d59'15.00"E, 57d28'23.20"N)
Center      (  554075.000, 6498580.000) ( 24d55'50.87"E, 58d37'25.43"N)
Band 1 Block=37013x1 Type=Float32, ColorInterp=Gray

And on the (aborted) output file:

Driver: GTiff/GeoTIFF
Files: ndvitest_commandline.tif
Size is 37100, 25800
Coordinate System is:
PROJCRS["Estonian Coordinate System of 1997",
    BASEGEOGCRS["EST97",
        DATUM["Estonia 1997",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4180]],
    CONVERSION["Estonian National Grid",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",57.5175539305556,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",24,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",59.3333333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",58,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",6375000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["northing (X)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (Y)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Topographic mapping (large scale)."],
        AREA["Estonia - onshore and offshore."],
        BBOX[57.52,20.37,60,28.2]],
    ID["EPSG",3301]]
Data axis to CRS axis mapping: 2,1
Origin = (369000.000000000000000,6635000.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
  PREDICTOR=2
Corner Coordinates:
Upper Left  (  369000.000, 6635000.000) ( 21d39'50.66"E, 59d49'52.42"N)
Lower Left  (  369000.000, 6377000.000) ( 21d48'50.21"E, 57d30'58.84"N)
Upper Right (  740000.000, 6635000.000) ( 28d16'31.78"E, 59d46'58.91"N)
Lower Right (  740000.000, 6377000.000) ( 28d 0' 5.93"E, 57d28'16.39"N)
Center      (  554500.000, 6506000.000) ( 24d56'23.65"E, 58d41'25.04"N)
Band 1 Block=37100x1 Type=Float32, ColorInterp=Gray
  NoData Value=-9999

Anyone have any clue what's going on here?

Improve this question
4
  • 1
    Please add the command with all parameters as text so it can be copied and pasted to the command line.
    – user30184
    Commented Mar 19, 2023 at 14:38
  • @user30184 Done
    – Clay
    Commented Mar 19, 2023 at 15:14
  • 1
    gdalwarp to VRT then gdal_translate to final tif gis.stackexchange.com/a/89549/2856
    – user2856
    Commented Mar 20, 2023 at 8:40
  • @user2856 it works! Many thanks to you. Post it as an answer so I can accept! Unless that's against the rules or something since it's linked from another answer.
    – Clay
    Commented Mar 20, 2023 at 9:21

2 Answers 2

Reset to default
5

Without your source data the best test I could make was to create as similar input file as you have with the gdal_create utility https://gdal.org/programs/gdal_create.html. I used GDAL 3.7.0dev on Windows

gdal_create -of GTiff -ot Float32 -outsize 37013 24284 -bands 1 -co blockxsize=37013 -co blockysize=1 -a_srs epsg:3301 -a_ullr 369010.000 6620000.000 739140.000 6377160.000 estonian_summer_2022_source.tif

gdalwarp -dstnodata -9999.0 -tr 10.0 10.0 -r med -te 369000.0 6377000.0 740000.0 6635000.0 -multi -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 estonian_summer_2022_source.tif estonian_summer_2022_target.tif

Creating output file that is 37100P x 25800L.
Processing estonian_summer_2022_source.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

So the command works for me with GDAL 3.7.0 with a similar, but empty, tiff file that you have.

Improve this answer
1
  • I didn't know about gdal_create, thanks for the tip. This works for me in Linux too, so somehow the problem is somewhere between the file itself and my Linux GDAL installation. (I also tried the real file in OSGeo shell in Windows and that worked.) I noticed that as command runs in Linux, the new file in the browser is continually growing in size to around 4 gigs, when it errors out. Whereas thru Windows OSGeo/QGIS, the new file stays at 1kb in size until it's done being processed. Something is different in the implementation here I think.
    – Clay
    Commented Mar 20, 2023 at 8:08
3

It's a well known and longstanding issue that gdalwarp doesn't deal with compression well. The solution is to gdalwarp without compression then gdal_translate with compression.

To avoid two lengthy processes, gdalwarp to VRT first, it's really quick, then gdal_translate with the -co compress option.

i.e.

$ gdalwarp etc... -of vrt input_file.tif output_file.vrt
$ gdal_translate -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 output_file.vrt output_file.tif

If using GDAL >=2 you can combine this into a single operation by writing the VRT to /vsistdout and piping that to gdal_translate and specifying /vsistdin as the input. For example:

gdalwarp -q  etc... -of vrt input_file.tif /vsistdout/ | gdal_translate -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9   /vsistdin/ output_file.tif
Improve this answer

Not the answer you're looking for? Browse other questions tagged or ask your own question.