How can I create a consistent <srs> element in a QGIS project's layer?

Question

From a standalone python script (i.e. not running within QGIS) I am manipulating QGIS a project file (*.qgs) to have its layers point to different geometry data.

That geometry data can be in different spatial reference systems (SRS), so I'd like to also adapt the <maplayer>'s <srs> child element, so that QGIS will know it and display it correctly in the layer properties dialog.

Unfortunately, (like many other things in the QGIS project files) this element seems to be unnecessarily redundant and QGIS seems to actually depend on redundant information for the display. The entry might look like this:

<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
<qgis projectname="" version="2.18.10">
  <!-- ... -->
  <projectlayers>
    <maplayer <!-- ... --> >
      <!-- ... -->
      <srs>
        <spatialrefsys>
          <proj4>+proj=longlat +ellps=intl +towgs84=72.438,345.918,79.486,1.6045,0.8823,0.5565,1.3746 +no_defs</proj4>
          <srsid>3482</srsid>
          <srid>4629</srid>
          <authid>EPSG:4629</authid>
          <description>Tahaa 54</description>
          <projectionacronym>longlat</projectionacronym>
          <ellipsoidacronym>intl</ellipsoidacronym>
          <geographicflag>true</geographicflag>
        </spatialrefsys>
      </srs>
      <!-- ... -->
    </maplayer>
  </projectlayers>
  <!-- ... -->
</qgis>

Having just the PostGIS/EPSG SRID, how do I get the rest of the information needed to create such an entry consistently? The SRS can differ but is known (by its SRID) at runtime.

I tried giving only the PostGIS/EPSG SRID by shortening the entry to

<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
<qgis projectname="" version="2.18.10">
  <!-- ... -->
  <projectlayers>
    <maplayer <!-- ... --> >
      <!-- ... -->
      <srs>
        <spatialrefsys>
          <srid>4629</srid>
        </spatialrefsys>
      </srs>
      <!-- ... -->
    </maplayer>
  </projectlayers>
  <!-- ... -->
</qgis>

but then QGIS will not display it correctly in the layer properties: expected:

(Also, while the data is being displayed, I'm not sure whether it would be projected correctly with that partial information.)

Where does QGIS take this information from when creating a new layer?

Answer 1

Do it like QGIS does it internally.

1. Create a QgsCoordinateReferenceSystem
2. Initialize with createFromId
3. use the objects writeXml method

Answer 2

QGIS reads CRS data from a SQLite db, on my machine this is in /usr/share/qgis/resources/srs.db

A rummage in sqlite3 shows several tables and a view (usual caveat, be careful with this!)

These tables (and view) seem to have most of the information shown in the XML snippet as stored in the .qgs file

In sqlite3 , .schema shows the tables. One of interest is tbl_srs, which has the following definition:

-- Don't run this SQL code.
-- It's just the table definition and would fail anyway.

CREATE TABLE tbl_srs (
  srs_id INTEGER PRIMARY KEY,
  description text NOT NULL,
  projection_acronym text NOT NULL,
  ellipsoid_acronym NOT NULL,
  parameters text NOT NULL,
  srid integer NOT NULL,
  auth_name varchar NOT NULL,
  auth_id varchar NOT NULL,
  is_geo integer NOT NULL
, deprecated boolean, noupdate boolean);

You can see the mapping from srs_id to auth_id

sqlite> .headers on
sqlite> .mode column
sqlite> select srs_id, description, auth_id from tbl_srs where srs_id=3482;
srs_id      description  auth_id   
----------  -----------  ----------
3482        Tahaa 54     4629      

sqlite> select srs_id, description, auth_id from vw_srs where auth_id=3482;
srs_id      description                     auth_id   
----------  ------------------------------  ----------
1440        NAD83(NSRS2007) / Arizona West  3482      

If you look at the QGIS CRS dialogs, auth_id is what is used. EPSG:3482 is something entirely different (NAD83(NSRS2007) / Arizona West)

You can also query the proj4 string. I'm pretty sure there's a way in pyqgis to extract individual parts from a proj4 string/CRS, or to convert it into WKT format. But a lot (if not all) of the fields you see in the project file XML are available from this database.

I believe that everything that QGIS knows about a given coordinate system is gathered from those tables / views.

SQLite support is a standard Python feature since Python 2.5 so it should be reasonably easy to implement.

Answer 3

I would guess all that information is available somewhere inside QGIS but another source would be to make a call (or two) to http://epsg.io, in your example you can get the proj4 string by parsing http://epsg.io/4629 to get

+proj=longlat +ellps=intl +towgs84=72.438,345.918,79.486,1.6045,0.8823,0.5565,1.3746 +no_defs 

and if you look at the response from http://epsg.io/4629.wkt

GEOGCS["Tahaa 54",
DATUM["Tahaa_54",
    SPHEROID["International 1924",6378388,297,
        AUTHORITY["EPSG","7022"]],
    TOWGS84[72.438,345.918,79.486,1.6045,0.8823,0.5565,1.3746],
    AUTHORITY["EPSG","6629"]],
PRIMEM["Greenwich",0,
    AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
    AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4629"]]

You can probably pick up the rest of the parameters.