WarpedVRT and SRC_METHOD GEOLOC_ARRAY

[PacoDu]

Hi !

I'm trying to use rasterio to reproject a netCDF file with lon lat grid coordinates to a COG file, but I can't find the appropriate way to call rasterio to reproduce the equivalent gdalwarp call:

gdalwarp -of COG -geoloc input.nc output.tiff 

As far as I understand I need to pass the option SRC_METHOD=GEOLOC_ARRAY to the GDALCreateGenImgProjTransformer2 to reproduce the equivalent of -geoloc option:
https://gdal.org/api/gdal_alg.html#_CPPv432GDALCreateGenImgProjTransformer212GDALDatasetH12GDALDatasetHPPc

But using the following code is not working (inspired by this test https://github.com/rasterio/rasterio/blob/master/tests/test_warpedvrt.py#L205), rasterio doesn't properly interpret the geolocation from the source netcdf and the output tiff is not reprojected:

    config = dict(
        GDAL_NUM_THREADS="ALL_CPUS",
    )
    transformer_options = { "SRC_METHOD": "GEOLOC_ARRAY" }

    with rasterio.Env(**config):
        with rasterio.open(src_path) as src:
            with WarpedVRT(
                src,
                # src_crs=CRS.from_epsg(4326),
                crs=CRS.from_epsg(3857),
                **transformer_options
            ) as vrt:
                output_profile = cog_profiles.get("raw")

                cog_translate(
                    vrt,
                    dst_path,
                    output_profile,
                    config=config,
                    in_memory=False,
                    quiet=True,
                    **options,
                )

If I set src_crs=CRS.from_epsg(4326) in WarpedVRT I get the following error:

rasterio._err.CPLE_AppDefinedError: Too many points (5625 out of 5625) failed to transform, unable to compute output bounds.

Here is the expected VRT Transformer output (generated with gdalwarp -of VRT -geoloc input.nc output.vrt):

<GDALWarpOptions>
(...)
    <Transformer>
      <ApproxTransformer>
        <MaxError>0.125</MaxError>
        <BaseTransformer>
          <GenImgProjTransformer>
            <SrcGeoLocTransformer>
              <GeoLocTransformer>
                <Reversed>0</Reversed>
                <Metadata>
                  <MDI key="LINE_OFFSET">0</MDI>
                  <MDI key="LINE_STEP">1</MDI>
                  <MDI key="PIXEL_OFFSET">0</MDI>
                  <MDI key="PIXEL_STEP">1</MDI>
                  <MDI key="SRS">GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]</MDI>
                  <MDI key="X_BAND">1</MDI>
                  <MDI key="X_DATASET">NETCDF:"data/input.nc":lon</MDI>
                  <MDI key="Y_BAND">1</MDI>
                  <MDI key="Y_DATASET">NETCDF:"data/input.nc":lat</MDI>
                </Metadata>
              </GeoLocTransformer>
            </SrcGeoLocTransformer>
            <DstGeoTransform>-80.7027849579179559,0.0425886528139265627,0,73.6016328141495961,0,-0.0425886528139265627</DstGeoTransform>
            <DstInvGeoTransform>1894.93631814359742,23.4804327896701679,0,1728.19819250262208,0,-23.4804327896701679</DstInvGeoTransform>
          </GenImgProjTransformer>
        </BaseTransformer>
      </ApproxTransformer>
    </Transformer>
(...)
</GDALWarpOptions>

It looks like I'm doing something wrong or maybe this is not supported by rasterio ?
Thanks for any feedback ! 🙏

GDAL version: GDAL 3.3.2, released 2021/09/01
rasterio version: 1.2.10
rio-cogeo version: 3.2.0

[vincentsarago]

@PacoDu would you be able to share the file, it would help I think

[PacoDu]

Sure here is the netcdf file download

[PacoDu]

Note that I'm facing this issue with multiple netcdf files that follow the CF-1.6 convention

[sgillies]

@PacoDu I'm sorry to be so late to reply. Support for transformer options like SRC_METHOD was added in Rasterio 1.3.0: #2195. You would need to upgrade.