Mention importance of PROJ grids and PROJ_NETWORK

[dugalh]

Expected behavior and actual behavior.

Hello. rasterio.warp.reproject() with apply_vertical_shift=True is not applying vertical shifts. When both src_crs and dst_crs have vertical components, I would expect a vertical shift to be applied. I see the same issue with rasterio.warp.transform().

FYI this behaviour does not seem to be consistent across different rasterio packages - the conda-forge packages work ok.

Steps to reproduce the problem.

import numpy as np
import rasterio as rio
from rasterio.warp import reproject

src_crs = 'EPSG:32634+5773'
dst_crs_list = ['EPSG:32634', 'EPSG:32634+4326', 'EPSG:32634+3855']

src_array = np.ones((1, 1, 1)) * 10
src_transform = rio.transform.from_origin(10, 10, 1, 1)

for dst_crs in dst_crs_list:
    dst_array, dst_transform = reproject(
        src_array, src_transform=src_transform, src_crs=src_crs, dst_crs=dst_crs,
        apply_vertical_shift=True
    )
    print(f'{src_crs}: {src_array.item(0):.4f}, {dst_crs}: {dst_array.item(0):.4f}')

This gives me::

EPSG:32634+5773: 10.0000, EPSG:32634: 10.0000
EPSG:32634+5773: 10.0000, EPSG:32634+4326: 10.0000
EPSG:32634+5773: 10.0000, EPSG:32634+3855: 10.0000

But I think it should give::

EPSG:32634+5773: 10.0000, EPSG:32634: 10.0000
EPSG:32634+5773: 10.0000, EPSG:32634+4326: 4.3727
EPSG:32634+5773: 10.0000, EPSG:32634+3855: 9.8771

which is what I get with a conda-forge rasterio package.

Environment Information

rasterio info:
  rasterio: 1.3.8
      GDAL: 3.6.4
      PROJ: 9.0.1
      GEOS: 3.11.1
 PROJ DATA: /home/dugalh/python-venv/py311-geospatial/lib/python3.11/site-packages/rasterio/proj_data
 GDAL DATA: /home/dugalh/python-venv/py311-geospatial/lib/python3.11/site-packages/rasterio/gdal_data

System:
    python: 3.11.0rc1 (main, Aug 12 2022, 10:02:14) [GCC 11.2.0]
executable: /home/dugalh/python-venv/py311-geospatial/bin/python3.11
   machine: Linux-5.15.90.1-microsoft-standard-WSL2-x86_64-with-glibc2.35

Python deps:
    affine: 2.4.0
     attrs: 23.1.0
   certifi: 2023.07.22
     click: 8.1.7
     cligj: 0.7.2
    cython: None
     numpy: 1.26.0
    snuggs: 1.4.7
click-plugins: None
setuptools: 59.6.0

Installation Method

rasterio-1.3.8-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl installed from PyPI with pip 22.0.2,

[sgillies]

@dugalh can you check the GDAL change log and see if there was a fix after 3.6.4?

[dugalh]

@sgillies I've looked into this some more - it is not a bug but a PROJ behaviour I was unaware of. The transformation grids needed for the vertical shift are not included the wheel, but can be downloaded by hand, or downloaded and cached on demand with the right settings. There are sections in the pyproj and PROJ docs describing this.

It does the vertical shift if I set the PROJ_NETWORK environment variable to 'ON' before calling warp.reproject or warp.transform:

import os
import numpy as np
import rasterio as rio
from rasterio.warp import reproject

os.environ.update(PROJ_NETWORK='ON')  # <--

src_crs = 'EPSG:32634+5773'
dst_crs_list = ['EPSG:32634', 'EPSG:32634+4326', 'EPSG:32634+3855']

src_array = np.ones((1, 1, 1)) * 10
src_transform = rio.transform.from_origin(10, 10, 1, 1)

for dst_crs in dst_crs_list:
    dst_array, dst_transform = reproject(
        src_array, src_transform=src_transform, src_crs=src_crs, dst_crs=dst_crs,
        apply_vertical_shift=True
    )
    print(f'{src_crs}: {src_array.item(0):.4f}, {dst_crs}: {dst_array.item(0):.4f}')

Which outputs:

EPSG:32634+5773: 10.0000, EPSG:32634: 4.3727
EPSG:32634+5773: 10.0000, EPSG:32634+4326: 4.3727
EPSG:32634+5773: 10.0000, EPSG:32634+3855: 9.8771

In the first line it seems the dst_crs of EPSG:32634 is being promoted to EPSG:32634+4326, which is different to the conda-forge rasterio behaviour:

EPSG:32634+5773: 10.0000, EPSG:32634: 10.0000
EPSG:32634+5773: 10.0000, EPSG:32634+4326: 4.3727
EPSG:32634+5773: 10.0000, EPSG:32634+3855: 9.8771

Replicating with pyproj, I get the first output with pyproj version 3.3.1 (PROJ version 8.2.0), and the second output with pyproj version 3.6.1 (PROJ version 9.3.0). So it seems the anomaly relates to the PROJ version, if that all makes sense...

[dugalh]

Perhaps there could be a note with the warp.reproject and warp.transform docs to make this more visible to users?