6

I'm looking for a method to shift the poles vertically in Mollweide projection using QGIS 3.

I have come across this website and I want to do something similar.

This is the default view of the Mollweide projection by Jason Davies. Default Mollweide projection

This is my take on Mollweide projection in QGIS 3 with minor alterations. enter image description here

This is Mollweide projection by Jason Davies with the poles shifted vertically enter image description here

How can I achieve the same result in QGIS 3?

3
  • 1
    Do you know where is it centered? I have created a new datum in proj, but seems to me that a pipeline in ogr2ogr must be possible. Also, is not easy to understand where to cut the Antarctica polygon: i.stack.imgur.com/HC92N.png Jan 10 '20 at 4:46
  • Thanks Gabriel. I want to center my map somewhere around these coordinates: 59.9375, 30.30861. Can you recommend some sources where I can learn about creating new datum in proj?
    – iamcyrk
    Jan 10 '20 at 10:48
  • I have just added an answer. If you have any doubt you can ask in its comments, or you if you are stuck at any stage of the process you can start a new question referring to this one. Jan 11 '20 at 4:01
3

Since Mollweide does not accept a latitude shift, what we can do is transform the coordinates of the geometries to a new geographic CRS with the Poles shift.

For that, I will register a custom geographic CRS in the PROJ database, and associate to it a Helmert Transformation.

PROJ registers all data in a SQLite database. I will edit it with DB Brower for SQLite.

Since I will do this on an OSGeo4W install over a Windows 64bit machine, I will open, in DB Browser, the C:\OSGeo4W64\share\proj\proj.db file.


First, I will create a new spherical datum. That can be done in the Browse Data tab, creating a new row in the geodetic_datum table:

1

Or it can be done in the Execute SQL tab, running the following query:

INSERT INTO "main"."geodetic_datum"
("auth_name", "code", "name", "ellipsoid_auth_name", "ellipsoid_code", "prime_meridian_auth_name", "prime_meridian_code", "area_of_use_auth_name", "area_of_use_code", "deprecated")
VALUES ('CUSTOM', '200001', 'A_Sphere', 'EPSG', '7048', 'EPSG', '8901', 'EPSG', '1262', 0);

2

Where:

  • CUSTOM:200001 will be the auth:code of our datum.

  • A_Sphere will be its name.

  • EPSG:7048 is the auth:code for the GRS 1980 Authalic Sphere ellipsoid.
    It is a sphere which surface is the same as the GRS 80 ellipsoid (R=6371007.0 m).

  • EPSG:8901 is the auth:code for Greenwich prime meridian.

  • EPSG:1262 is the auth:code for a World area of use.

And I will register its allias name with the following query:

INSERT INTO "main"."alias_name"
("table_name", "auth_name", "code", "alt_name", "source")
VALUES ('geodetic_datum', 'CUSTOM', '200001', 'A_Sphere', 'CUSTOM');

Then, I will create a new geographic coordinate system, with the following query:

INSERT INTO "main"."geodetic_crs"
("auth_name", "code", "name", "type", "coordinate_system_auth_name", "coordinate_system_code", "datum_auth_name", "datum_code", "area_of_use_auth_name", "area_of_use_code", "deprecated")
VALUES ('CUSTOM', '200101', 'Shifted_Poles', 'geographic 2D', 'EPSG', '6422', 'CUSTOM', '200001', 'EPSG', '1262', 0);

Where:

  • CUSTOM:200101 will be the auth:code of our new geodetic system.

  • Shifted_Poles will be its name.

  • EPSG:6422 is the auth:code for ellipsoidal 2D coordinate systems (our sphere is a particular case of an ellipsoid).

  • CUSTOM:200001 will be the daum associated to our coordinate system.

  • 0 means for it is not deprectaed.


Finally, I will register a new Helmert transformation from our coordinate system to WGS 84. In this transformation we will ask to rotate 30.30861 degrees around the Y axis, and -59.9375 degrees around the Z axis.
Z axis in the ellipsoid is coincident with the minor axis (line of Poles). X axis belongs in the equator plane and points to the prime meridian. And Y axis belongs in the equator plane and points to longitude 90 degrees. Rotations are clockwise.

INSERT INTO "main"."helmert_transformation_table"
("auth_name", "code", "name", "method_auth_name", "method_code", "source_crs_auth_name", "source_crs_code", "target_crs_auth_name", "target_crs_code", "area_of_use_auth_name", "area_of_use_code", "accuracy", "tx", "ty", "tz", "translation_uom_auth_name", "translation_uom_code", "rx", "ry", "rz", "rotation_uom_auth_name", "rotation_uom_code", "scale_difference", "scale_difference_uom_auth_name", "scale_difference_uom_code", "operation_version", "deprecated")
VALUES ('CUSTOM', '200201', 'Poles_Shift', 'EPSG', '9607', 'CUSTOM', '200101', 'EPSG', '4326', 'EPSG', '1262', 0.0, 0.0, 0.0, 0.0, 'EPSG', '9001', 0.0, 30.30861, -59.9375, 'EPSG', '9102', 0.0, 'EPSG', '9202', '1.0.0', 0);  

Where:

  • CUSTOM:200201 will be the auth:code of our transformation.

  • Poles_Shift will be its name.

  • EPSG:9607 is the auth:code for a Coordinate Frame Rotation (geog2D domain).

  • CUSTOM:200101 is the source CRS auth:code.

  • EPSG:4326 is the target CRS auth:code.

  • The first 0.0 parameter is the accuracy (not true value in our case, just informative).

  • 0.0, 0.0, 0.0, 'EPSG, '9001' are the translation parameters, in the meters auth:code.

  • 0.0, 30.30861, -59.9375, 'EPSG', '9102' are the rotation parameters, in the degrees auth:code.

  • 0.0, 'EPSG', '9202' are the scale parameter, in p.p.m. auth:code.

  • 1.0.0 is our version of this transformation.


Ok, don't forget to Write Changes in the database! That was our main work. We have defined a new spherical datum, a new geographic coordinate system based on that datum, and a new transformation between our system and WGS84 system.

What that transformation means for: The latitude 0 degrees, longitude 0 degrees of our system, will be approximately in the latitude 30.30861 degrees, longitude 59.9375 degrees in the WGS84 system (not excactly because the ellipsoidal shape of its datum).

Let's see how projinfo program interpret our system. Run the C:\OSGeo4W64\OSGeo4W.bat file. In the console opened:

C:\>projinfo CUSTOM:200101
PROJ.4 string:
+proj=longlat +R=6371007 +towgs84=0,0,0,0,-109110.996,215775,0 +no_defs +type=crs

WKT2_2018 string:
GEOGCRS["Shifted_Poles",
    DATUM["A_Sphere",
        ELLIPSOID["GRS 1980 Authalic Sphere",6371007,0,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]],
    ID["CUSTOM",200101]]

Looks ok to me. +towgs84 parameters are converted in seconds there.

Let's see how the transformation from our system to WGS84 must work:

C:\>projinfo -s CUSTOM:200101 -t EPSG:4326
Candidate operations found: 1
-------------------------------------
Operation n┬░1:

CUSTOM:200201, Poles_Shift, 0.0 m, World

PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=longlat +R=6371007 +step +proj=push +v_3 +step +proj=cart +R=6371007 +step +proj=helmert +x=0 +y=0 +z=0 +rx=0 +ry=109110.996 +rz=-215775 +s=0 +convention=coordinate_frame +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

WKT2_2018 string:
COORDINATEOPERATION["Poles_Shift",
    VERSION["1.0.0"],
    SOURCECRS[
        GEOGCRS["Shifted_Poles",
            DATUM["A_Sphere",
                ELLIPSOID["GRS 1980 Authalic Sphere",6371007,0,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["CUSTOM",200101]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    METHOD["Coordinate Frame rotation (geog2D domain)",
        ID["EPSG",9607]],
    PARAMETER["X-axis translation",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8605]],
    PARAMETER["Y-axis translation",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8606]],
    PARAMETER["Z-axis translation",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8607]],
    PARAMETER["X-axis rotation",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8608]],
    PARAMETER["Y-axis rotation",30.30861,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8609]],
    PARAMETER["Z-axis rotation",-59.9375,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8610]],
    PARAMETER["Scale difference",0,
        SCALEUNIT["parts per million",1E-06],
        ID["EPSG",8611]],
    OPERATIONACCURACY[0.0],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]],
    ID["CUSTOM",200201]]

That's nice! We have now PROJ.4 string and WKT2_2018 definitions for our system, and PROJ pipeline and WKT2_2018 definition for our transformation. However, now that we know them, we no longer need them... We can refer to our system as CUSTOM:200101 from now.


Let's transform our GeoPackages:

C:\GA\GIS\Pruebas>ogr2ogr -t_srs CUSTOM:200101 continentes2.gpkg continentes.gpkg

C:\GA\GIS\Pruebas>ogr2ogr -t_srs CUSTOM:200101 grilla2.gpkg grilla.gpkg

And load them in QGIS:

3

Note that I am not reprojecting the canvas (Project -> Properties -> CRS -> No projection). We do not want to reproject the geometries again. Doing it will show them approximately where they must be in the Earth.

What I will do is Set the EPSG:4326 CRS to the new layers.

In that way, we are saying to QGIS, interpret those geometries as if them were in WGS84.

I'm sorry about our continents that I must delete, because they crosses the timeline. What I would do, is create a world bounding box in this system and transform it to WGS84 to see where our new timeline crosses them, and clip the geometries there.


Once the layers CRS is set, reproject them to Mollwide:

4


External References:

5
  • About how to handle the geometries, also with the use of the pipeline without registering the system in PROJ database, I have asked a new question: gis.stackexchange.com/questions/347157/… Jan 12 '20 at 0:09
  • Thanks a lot, Gabriel. I have followed your instructions up to this line: "We can refer to our system as CUSTOM:200101 from now." Everything works like a charm. Now I have a bit of a trouble with this line: "Let's transform our GeoPackages: C:\GA\GIS\Pruebas>ogr2ogr -t_srs CUSTOM:200101 continentes2.gpkg continentes.gpkg C:\GA\GIS\Pruebas>ogr2ogr -t_srs CUSTOM:200101 grilla2.gpkg grilla.gpkg" How do I transform GeoPackages? Why is there a reference to ogr2ogr library? What is the general purpose of these 2 strings of code?
    – iamcyrk
    Jan 16 '20 at 12:18
  • 1
    @iamcyrk, Hi, you are welcome. We need ogr2ogr GDAL program to transform the coordinates of the geometries. PROJ.7 will be released soon, and QGIS 3.12 will adopt custom pipelines internally. But for now we must transform the coordinates with ogr2ogr. It is included with all QGIS installs. For now, maybe a Custom CRS with the PROJ.4 string can be created in QGIS, and use that CRS to export the layers from QGIS. Since we have found PROJ.4 string and PROJ.6 pipeline for the transformed CRS, we have more ways to do that. Check my question post referred in my previous comment. Jan 16 '20 at 12:41
  • Thank you. But I still have some issues. I managed to create continentes2.gpkg file. But when I try to run this string C:\OSGeo4W64\bin>ogr2ogr -t_srs CUSTOM:200101 continentes2.gpkg continentes.gpkg I receive a failure notice: Unable to open datasource 'continentes.gpkg' with the following drivers. -> JP2ECW' -> OCI' -> SOSI' -> PCIDSK' and the list goes on up to -> `HTTP' What can be the cause and what should I do now?
    – iamcyrk
    Jan 21 '20 at 13:36
  • @iamcyrk, In my answer there is a link to the ogr2ogr manual page. continentes2.gpkg is the destination file. The source file is continentes.gpkg. Jan 21 '20 at 14:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.