Cannot extract by mask using gdalwarp

Question

I am trying to extract data from a raster file using a Shapefile mask. QGIS suggests the following command:

gdalwarp -of GTiff -cutline SanFernando_basin.shp -cl SanFernando_basin -crop_to_cutline SFBimg.tif SFBout.tif

However, when I run that command in terminal, get the following error message:

ERROR 1: PROJ: proj_create_operations: SQLite error on SELECT source_crs_auth_name, source_crs_code, target_crs_auth_name, target_crs_code, cov.auth_name, cov.code, cov.table_name, area.south_lat, area.west_lon, area.north_lat, area.east_lon, ss.replacement_auth_name, ss.replacement_code FROM coordinate_operation_view cov JOIN area ON cov.area_of_use_auth_name = area.auth_name AND cov.area_of_use_code = area.code LEFT JOIN supersession ss ON ss.superseded_table_name = cov.table_name AND ss.superseded_auth_name = cov.auth_name AND ss.superseded_code = cov.code AND ss.superseded_table_name = ss.replacement_table_name AND ss.same_source_target_crs = 1 WHERE ((source_crs_auth_name = ? AND source_crs_code = ? AND target_crs_auth_name = ? AND target_crs_code = ?) OR (source_crs_auth_name = ? AND source_crs_code = ? AND target_crs_auth_name = ? AND target_crs_code = ?)) AND cov.deprecated = 0 AND cov.auth_name = ? ORDER BY pseudo_area_from_swne(south_lat, west_lon, north_lat, east_lon) DESC, (CASE WHEN accuracy is NULL THEN 1 ELSE 0 END), accuracy: no such column: ss.same_source_target_crs
ERROR 6: Cannot find coordinate operations from `EPSG:4326' to `GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]'
ERROR 1: PROJ: proj_create_operations: SQLite error on SELECT source_crs_auth_name, source_crs_code, target_crs_auth_name, target_crs_code, cov.auth_name, cov.code, cov.table_name, area.south_lat, area.west_lon, area.north_lat, area.east_lon, ss.replacement_auth_name, ss.replacement_code FROM coordinate_operation_view cov JOIN area ON cov.area_of_use_auth_name = area.auth_name AND cov.area_of_use_code = area.code LEFT JOIN supersession ss ON ss.superseded_table_name = cov.table_name AND ss.superseded_auth_name = cov.auth_name AND ss.superseded_code = cov.code AND ss.superseded_table_name = ss.replacement_table_name AND ss.same_source_target_crs = 1 WHERE ((source_crs_auth_name = ? AND source_crs_code = ? AND target_crs_auth_name = ? AND target_crs_code = ?) OR (source_crs_auth_name = ? AND source_crs_code = ? AND target_crs_auth_name = ? AND target_crs_code = ?)) AND cov.deprecated = 0 AND cov.auth_name = ? ORDER BY pseudo_area_from_swne(south_lat, west_lon, north_lat, east_lon) DESC, (CASE WHEN accuracy is NULL THEN 1 ELSE 0 END), accuracy: no such column: ss.same_source_target_crs
ERROR 6: Cannot find coordinate operations from `EPSG:4326' to `EPSG:4326'

I do not know how to debug this error. Both the shapefile and the geotiff are in the same coordinate system (WGS84 / EPSG:4326), but I am not trying to reshape them. Just to apply the cutline to extract the relevant data.

Answer 1

I think the problem is with the way you have specified your OGR datasource i.e. ShapeFile. Test with:

gdalwarp -of GTiff -cutline /path/to/shapefile-directory -cl SanFernando_basin -crop_to_cutline SFBimg.tif SFBout.tif

Answer 2

By looking at the Proj source code the column name "same_source_target_crs" appears in three places: https://github.com/OSGeo/PROJ/search?q=same_source_target_crs

The new column definition has been added 2020-05-25 in commit https://github.com/OSGeo/PROJ/commit/4d9fa77bef599353d3809541cb7e6adb8ec12525

The commit has a description

Database: add a same_source_target_crs column to supersession table

This is in preparation for EPSG 9.8.11 import that supersedes a number of grid transformation for US transformations, but the superseded and replacement transformations don't operate on the same (source_crs, target_crs), which is a bit weird. So in that situation, ignores the supersession.

Obviously your QGIS has a Proj version that contains the above mentioned commit, but a proj.db database that has been created by an older version of Proj. Proj.db is a SQLite database that is used for storing all the parameters that Proj needs to know.

It could be enough to find a newer proj.db file and replace the old one with that. A safer alternative could be to install the whole Proj again, including the contents of the PROJ_LIB directory. If the same issue happens again there may be a bug in the conda packaging.