[QGIS-Developer] Dealing with datum shift grids in QGIS

[Even Rouault]

Hi,

some updates:

> I guess grid_alternatives could receive additional columns, which would be
> used if package_name is NULL, perhaps something like:
> 
> open_license BOOLEAN
> download_url TEXT
> direct_download BOOLEAN
> 

Done

> 
> > One final thought - how about customization of proj.db by ordinary
> > users? In QGIS users are allowed to define their custom CRS based on
> > Proj string - is something like that going to be supported?
> 
> The issue is that a PROJ string is a rather low expressive way of defining a
> CRS, but yes that something I've in a corner of my mind, to allow defining
> custom CRS entries by PROJ strings or WKT strings, without going to the
> full relational structure.

Done.

The geodetic_crs and projected_crs tables have now a 'text_definition' where users can put PROJ strings or WKT strings

> 
> > E.g. to
> > have system proj.db and user proj.db with customizations?
> 

Done (for the reading part by PROJ)

There's now a

int PROJ_DLL proj_context_set_database_path(PJ_CONTEXT *ctx,
                                            const char *dbPath,
                                            const char *const *auxDbPaths)

where:
- dbPath is the path to the 'system' proj.db (or NULL to let it find it at its default location)
- auxDbPaths is a list of auxiliary databases. They can potentially
refer to objects of the main database (if you don't need standalone integrity)

Demo:

1) Create an empty DB (custom.db) with the table and view structure of the system DB (proj.db)

echo "SELECT sql || ';' FROM sqlite_master WHERE type IN ('table', 'view');" | sqlite3 data/proj.db | sqlite3 custom.db

2) Fill it with a custom CRS definition expressed as old-school PROJ.4 string

echo "INSERT INTO projected_crs (auth_name, code, name, coordinate_system_auth_name , coordinate_system_code , geodetic_crs_auth_name, geodetic_crs_code, conversion_auth_name, conversion_code, area_of_use_auth_name, area_of_use_code, text_definition , deprecated) VALUES ( 'QGIS', 'MY-CRS-CODE', 'My custom CRS', NULL, NULL, NULL, NULL, NULL, NULL,  NULL, NULL, '+proj=utm +zone=32 +ellps=GRS80 +units=us-ft +towgs84=1,2,3', 0);" | sqlite3 custom.db

3) Query the DB:

$ src/projinfo --main-db-path data/proj.db --aux-db-path custom.db QGIS:MY-CRS-CODE -o PROJ,PROJ4,WKT2_2015

PROJ string: 
+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=utm +zone=32 +ellps=GRS80 +step +proj=unitconvert +xy_in=m +z_in=m +xy_out=us-ft +z_out=us-ft

PROJ.4 string: 
+proj=utm +zone=32 +ellps=GRS80 +towgs84=1,2,3,0,0,0,0 +units=us-ft

WKT2_2015 string: 
PROJCRS["My custom CRS",
    BASEGEODCRS["unknown",
        DATUM["Unknown based on GRS80 ellipsoid",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]]],
    CONVERSION["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["US survey foot",0.304800609601219]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["US survey foot",0.304800609601219]],
    ID["QGIS","MY-CRS-CODE"]]


Pedantic note: the new PROJ string format doesn't capture the towgs84, which is a legacy way of embedding a coordinate transformation. Only the PROJ.4 output does that for the traditionnalists. However internally PROJ has built a BoundCRS object to model the towgs84 transformation that was in the text_definition, so for PROJ.4 output it can recover it, and for coordinate operations, it can use it as well. Not sure I'm very clear here :-), but in short, things should just work.

4) Test a reprojection

$ src/projinfo --main-db-path data/proj.db --aux-db-path custom.db -s QGIS:MY-CRS-CODE -t EPSG:32633 -o PROJ

PROJ string: 
+proj=pipeline +step +proj=unitconvert +xy_in=us-ft +z_in=us-ft +xy_out=m +z_out=m +step +inv +proj=utm +zone=32 +ellps=GRS80 +step +proj=cart +ellps=GRS80 +step +proj=helmert +x=1 +y=2 +z=3 +step +inv +proj=cart +ellps=WGS84 +step +proj=utm +zone=33 +ellps=WGS84


Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com