using `cog_translate` on netCDF4 with improper georeferencing

ryanjdillon commented 5 years ago

I am trying to use cog_translate to convert some netCDF4 datasets to COG, and it appears these files do not have proper georeferencing.

I have modified cogeo.py to accept vrt_params, so that I can manually pass the src_crs and src_transform via the WarpedVRTReaderBase mixin of WarpedVRT, as you can see here:

https://github.com/ryanjdillon/rio-cogeo/blob/cog_translate_vrt_params/rio_cogeo/cogeo.py#L87-L99

Then I'm passing the following vrt_params:

vrt_params = dict(src_transform=Affine(gdal_geotransform), src_crs='EPSG:4326')

This get's me past the CRS is None error, but it doesn't seem to acknowledge the Affine transform I pass to src_transform, resulting in the following error:

CPLE_AppDefinedError: The transformation is already "north up" or a transformation between pixel/line and georeferenced coordinates cannot be computed for /home/ryan/data/webstep/kvt/NVE_Sorlandet/NVE_Sorlandet_wslev0_xyt_synth.nc4. There is no affine transformation and no GCPs. Specify transformation option SRC_METHOD=NO_GEOTRANSFORM to bypass this check.

Any suggestions on how I might get this to work?

vincentsarago commented 5 years ago

@ryanjdillon can you try transform instead of src_transform ? https://github.com/cogeotiff/rio-cogeo/blob/webOpt/rio_cogeo/cogeo.py#L141

see: https://github.com/mapbox/rasterio/blob/a84d5e1b2b00999e7c54cd47abf5443063fd7683/rasterio/_warp.pyx#L733-L742

ryanjdillon commented 5 years ago

Excellent. I'll report back as soon as I've tried that out. Thanks!

Den fre. 5. apr. 2019, 16:18 skrev Vincent Sarago <notifications@github.com

:

@ryanjdillon https://github.com/ryanjdillon can you try transform instead of src_transform ? https://github.com/cogeotiff/rio-cogeo/blob/webOpt/rio_cogeo/cogeo.py#L141

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/cogeotiff/rio-cogeo/issues/74#issuecomment-480292397, or mute the thread https://github.com/notifications/unsubscribe-auth/AA39Gv8T72PKy_nfaR93XtF17TUwujv0ks5vd1slgaJpZM4ce6K0 .

vincentsarago commented 5 years ago

@ryanjdillon did it fix the problem ?

ryanjdillon commented 5 years ago

I was just about to report back. That didn't seem to fix things unfortunately.

It looks like dst_crs and dst_transform are depreciated in favor of crs and transform, so I'm passing things with my modification like this


    vrt_params = dict(
        src_transform=nc_affine, src_crs=crs,
        transform=nc_affine, crs=crs,
        )

    cogeo.cog_translate(
        src_path=band_fp,
        dst_path=dst_fp,
        dst_kwargs=config.GDAL_TRANSLATE_PROFILE,
        vrt_params=vrt_params,
        indexes=[i,],
        overview_resampling='nearest',
        config=config.GDAL_CONFIG
        )

Which yields the same error. If I pass width and height, as used in the Affine instance, I get another error:

ObjectNullError: Pointer 'hDS' is NULL in 'GDALSetProjection'.

And that seems to be related to the output file not existing. It sounded as though this could have been due to the slightly outdated Debian packages I'm using, so I compiled everything from source, which didn't do anything for me.

gdalinfo --version

GDAL 2.4.1, released 2019/03/15

nc-config --all

This netCDF 4.6.3 has been built with the following features:

--cc -> gcc --cflags -> -I/usr/local/include -I/home/ryan/opt/hdf5-1.10.5/hdf5/include -I/home/ryan/opt/hdf-4.2.14/hdf4/include --libs -> -L/usr/local/lib -lnetcdf

--has-c++ -> no --cxx ->

--has-c++4 -> no --cxx4 ->

--has-fortran -> no --has-dap -> yes --has-dap2 -> yes --has-dap4 -> yes --has-nc2 -> yes --has-nc4 -> yes --has-hdf5 -> yes --has-hdf4 -> yes --has-logging -> no --has-pnetcdf -> no --has-szlib -> no --has-cdf5 -> yes --has-parallel4 -> no --has-parallel -> no

--prefix -> /usr/local --includedir -> /usr/local/include --libdir -> /usr/local/lib --version -> netCDF 4.6.3

ryanjdillon commented 5 years ago

I've tried to use cog_translate() as I described above on a file that is correctly georeferenced (Copernicus global NDVI from their Land service), and that gives me the same error.

However, gdalwarp seems to work by doing the following:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:4326 -of GTIFF 'NETCDF:"c_gls_NDVI_201808110000_GLOBE_PROBAV_V2.2.1.nc":NDVI' ./test.tiff

Whereas my NetCDF file lacking georeferencing throws the same error with gdalwarp as with cog_translate.

vincentsarago commented 5 years ago

@ryanjdillon so if I understand correctly this might be a Rasterio Error and not specific with rio-cogeo.

Can you check if you are at least able to read the dataset with a simple rasterio code ?

ryanjdillon commented 5 years ago

It seems like it is a bit of a combo. I got it to work when I opened the subdataset with a rasterio formatted band path, rather than gdal (i.e. src_path='HDF5:/path/to/file.nc://NDVI' instead of src_path='NETCDF:"/path/to/file.nc":NDVI').

It was also necessary to pass vrt_params (including src_crs, src_transform, crs, transform, width, and height).

I also tried using cog_validate on the generated tiff, and it bugged out when src.get_tag_item("BLOCK_OFFSET_0_0", "TIFF", bidx=1) was None.

I could submit a PR if you want to look over / incorporate these changes.

vincentsarago commented 5 years ago

🤔 thanks @ryanjdillon,

I also tried using cog_validate on the generated tiff, and it bugged out when src.get_tag_item("BLOCK_OFFSET_0_0", "TIFF", bidx=1) was None.

Can you checkout what https://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/validate_cloud_optimized_geotiff.py gives you ?

I could submit a PR if you want to look over / incorporate these changes

Allowing vrt_params as input might complicated a lot the code (especially with the changes from https://github.com/cogeotiff/rio-cogeo/pull/62). That said, I think we could document it and tell how to generate a proper temporary file from the netCDF4 before passing it to rio-cogeo.

I'm no expert with netCDF so I'm may be missing some knowledge here

ryanjdillon commented 5 years ago

I'll take a look at that as soon as possible.

Once I have a clear handle on things I can contribute an example to the docs or something.

ryanjdillon commented 5 years ago

I managed to get things to work by first converting the NetCDF to GeoTIFF using GDAL, and then doing the COG conversion with rio_cogeo.

I first created a VRT file with gdal_translate, passing the corner coords as GCPs. Then I used gdalwarp to convert from NetCDF to GeoTIFF. Passing the same GCPs to rasterio did not seem to correctly preserve the georeference in the resulting GeoTIFF. I am guessing this has to do with the files not conforming to a NetCDF "CF" standard.

The resulting COG produced by rio_cogeo passes the validation test in the osgeo script you linked to.

vincentsarago commented 4 years ago

@ryanjdillon Long time no see 😄 I finally had to work with netcdf recently and thanks to some recent change in rio-cogeo we can now pass Memory or VRT files to cog_translate.

Latest version

def translate(
src_path,
cog_path,
cogeo_profile="deflate",
cogeo_options={"blockxsize": 128, "blockysize": 128, "predictor": 2, "zlevel": 7},
):
"""Convert a netcdf to COG."""
config = dict(GDAL_NUM_THREADS="ALL_CPUS", GDAL_TIFF_OVR_BLOCKSIZE="128")

with rasterio.Env(**config):
    with rasterio.open(src_path) as src_dst:
        with WarpedVRT(
            src_dst,
            crs="epsg:4326",
            warp_mem_limit=0.1,   # Reduce chunk size, see: https://github.com/OSGeo/gdal/issues/1989
            num_threads=8
        ) as vrt:
            profile = cog_profiles.get(cogeo_profile)
            profile.update(cogeo_options)
            cog_translate(
                vrt,
                cog_path,
                profile,
                config=config,
                in_memory=True,
                quiet=True,
            )

First version

def translate(
src_path,
cog_path,
cogeo_profile="deflate",
cogeo_options={"blockxsize": 128, "blockysize": 128, "predictor": 2, "zlevel": 7},
):
"""Convert a netcdf to COG."""
config = dict(GDAL_NUM_THREADS="ALL_CPUS", GDAL_TIFF_OVR_BLOCKSIZE="128")

with rasterio.Env(**config):
    with rasterio.open(src_path) as src_dst:
        dst_transform, dst_width, dst_height = calculate_default_transform(
            src_dst.crs, "epsg:4326", src_dst.width, src_dst.height, *src_dst.bounds
        )
        meta = dict(
            driver="GTiff",
            dtype=src_dst.dtypes[0],
            count=src_dst.count,
            height=dst_height,
            width=dst_width,
            crs="epsg:4326",
            transform=dst_transform,
            nodata=-999,
            tiled=True,
            compress="deflate",
            blockxsize=256,
            blockysize=256,
        )
        with MemoryFile() as memfile:
            with memfile.open(**meta) as mem:
                reproject(
                    rasterio.band(src_dst, 1),
                    rasterio.band(mem, 1),
                    resampling=rasterio.enums.Resampling.nearest,
                    warp_mem_limit=1,  # Reduce chunk size, see: https://github.com/OSGeo/gdal/issues/1989
                    num_threads=8,
                )
                mem.update_tags(**src_dst.tags(1))
                profile = cog_profiles.get(cogeo_profile)
                profile.update(cogeo_options)
                cog_translate(
                    mem,
                    cog_path,
                    profile,
                    config=config,
                    in_memory=True,
                    quiet=True,
                )

since we can use MemoryFile or VRT, we can now fix the data upfront and pass it to cog_translate so I think we can close this issue.

ryanjdillon commented 4 years ago

Excellent. Thanks for the ping.

cogeotiff / rio-cogeo

using `cog_translate` on netCDF4 with improper georeferencing #74