[FR] Outputing Cloud Optimized Geotiffs (COG)

[ConnectedSystems]

Looking at how Rasters interacts with GDAL, it seems that I should be able to output a COG.

Rasters.write(
            mask_path,
            a_raster;
            ext=".tiff",
            source="gdal",
            driver="COG"
        )

As I understand it, only CreateCopy() is supported by GDAL but this seems to be handled by Rasters:

https://github.com/ConnectedSystems/Rasters.jl/blob/e97ddaf5dfa8febbebd2d08fea2ee73673ff0307/ext/RastersArchGDALExt/gdal_source.jl#L382-L391

However, it appears that a true COG is not created, despite passing the validation based on this python script

The resulting geotiff is not pixel interleaved and does not have a defined OVR_RESAMPLING_ALG - maybe this has something to do with it, maybe not. Wondering if someone can give me any pointers.

But either way, producing COGs would be nice to be able to do with Rasters.jl

I've included the file created with Rasters.jl in case relevant slopes_api_output.zip

[rafaqz]

Its meant to work if you specify "COG", but yes it probably needs some more testing.

I think the key step here to have a test written in Julia of what you mean by "true COG". If that python script passes I'm unlikely to do any better.

If you can write that we can put it in the tests, and I can do the part of actually fixing the output.

[rafaqz]

Yes GDAL and so Rasters can only create these types directly:

https://github.com/ConnectedSystems/Rasters.jl/blob/e97ddaf5dfa8febbebd2d08fea2ee73673ff0307/ext/RastersArchGDALExt/gdal_source.jl#L8

And otherwise its via copy. I have no idea why GDAL is like that. But "COG" is one driver where we create a tiff first, then copy it to COG. It's quite possible I have missed some things in that process, GDAL is kind of arcane and as you say things seem to work fine on the surface with the current code.

[ConnectedSystems]

Its meant to work if you specify "COG", but yes it probably needs some more testing.

I think the key step here to have a test written in Julia of what you mean by "true COG". If that python script passes I'm unlikely to do any better.

If you can write that we can put it in the tests, and I can do the part of actually fixing the output.

Part of the issue is that I'm not able to discern why the file I produce with Rasters.jl is not correctly interpreted by ArcGIS and related web maps as a COG.

For example, this file, which is linked in this example page

I ran gdalinfo on it as well as running that validation script (found this streamlit app that wraps around it), and nothing stands out to me as to why ArcGIS Online treats this as a COG, but not the one produced by Rasters.jl

[rafaqz]

Ah right. Probably I'm also unlikely to find anything. For background I'm finishing a PhD and snowed under with everything else in these packages, so I'm unlikely to chase after something underdefined like this before next year.

But, if you can find what the problem is exactly I will try to fix it. You may need to dig into the GDAL docs and see if the Rasters create/copy implementation does everything it needs to.

Another approach is to do it all manually using GDAL or ArchGDAL until it works, and then we can work off that successful MWE.

[ConnectedSystems]

Yeah, I may have to resort to using something cogger as a backup.

Will dig into this issue when I can though.

All the best with your PhD!

[rafaqz]

Sure, may be pragmatic. But do try to do it directly with ArchGDAL, maybe it will just work and should only take half hour to write.

[ConnectedSystems]

Sure, may be pragmatic. But do try to do it directly with ArchGDAL, maybe it will just work and should only take half hour to write.

I have taken a look already and yes it should be possible, but I'd prefer to get something working first then make it seamless. Thanks for all your input!

[rafaqz]

Yea I'm talking from a personal need to see a GDAL example that works so I can close this issue! I'm sure just using cogger is better for getting your actual project done

[ConnectedSystems]

@rafaqz

This works for our purposes:

import ArchGDAL as AG

test_data = AG.read("path_to_some.tiff")

AG.write(test_data, "archgdal_cog.tif"; driver=AG.getdriver("COG"))

The issue seems to be something to do with the default compression method. LZW/DEFLATE seems to work but ZSTD does not (which happens to be what many COG creation processes use, including cogger).

I suspect ArcGIS is using an older GDAL version or otherwise following an older standard for COGs.

I've attached a set of data used to explore this with a README in case of interest. I had a quick look at some point but it would be nice to be able to pass options directly to GDAL with Rasters.write(). If I can already do that, then I must have missed it when reading the docs/code.

cog_create_issue (2).zip

[ConnectedSystems]

Just for extra clarity:

As far as I can tell, what Rasters/ArchGDAL creates is a COG, just not one that is supported by ArcGIS. If one needs a COG that can be used with ArcGIS, use a different compression method.

For anyone else who struggled with this, the below is what we ended up going with. LZW and DEFLATE compression both seem to work, with LZW being slightly faster for our purposes.

Rasters.write(
      mask_path,
      raster_data;
      ext=".tiff",
      source="gdal",
      driver="COG",
      options=Dict(
          "COMPRESS"=>"LZW",
          "SPARSE_OK"=>"TRUE",
      )
)