7

I’m trying to plot a shapefile of points that correspond for all impact craters on Pluto. I can't get the points to math with the impact craters from the global mosaic of Pluto's surface. In the best o QGIS has a CRS named GCS_Pluto_2000, but when I set the global mosaic, the shape of points and the project in QGIS to the same CRS, the points are plotted completely wrong.

I might be missing something but I can’t guess what.

The link for the global mosaic and the CSV with the points:

https://astrogeology.usgs.gov/search/map/Pluto/NewHorizons/Pluto_NewHorizons_Global_Mosaic_300m_Jul2017

https://astrogeology.usgs.gov/search/map/Pluto/Research/Craters/Craters_PlutoCharon_System_Robbins

3
  • 1
    Longitude values are in 0-360. If the image is +/-180, you'll have to convert the values into the same range. – mkennedy Mar 18 '20 at 18:00
  • I have already done it. But the points still do not match the craters . Some times it misses for a few kilometers – caio villaca Mar 18 '20 at 18:03
  • 1
    Hmmm, it could be that one is using the equivalent of "geographic" coordinates and the other is planetographic. If they're off up to a few km, it wouldn't be a 0-360 positive west vs east at least. You may need to dig into the metadata--see if they're both using the Pluto 2000 definition, etc. – mkennedy Mar 18 '20 at 18:10
11

The mosaic has a Projected (Equidistant Cyclindrical) CRS. We can check it with:

gdalinfo Pluto_NewHorizons_Global_Mosaic_300m_Jul2017_8bit.tif

Prints:

Driver: GTiff/GeoTIFF
Files: Pluto_NewHorizons_Global_Mosaic_300m_Jul2017_8bit.tif
Size is 24888, 12444
Coordinate System is:
PROJCRS["SimpleCylindrical Pluto",
    BASEGEOGCRS["GCS_Pluto",
        DATUM["D_Pluto",
            ELLIPSOID["Pluto",1188300,0,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Reference_Meridian",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["Equidistant Cylindrical",
        METHOD["Equidistant Cylindrical",
            ID["EPSG",1028]],
        PARAMETER["Latitude of 1st standard parallel",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Longitude of natural origin",180,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2
Origin = (-3733200.000000000000000,1866600.000000000000000)
Pixel Size = (300.000000000000000,-300.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-3733200.000, 1866600.000) (  0d 0' 7.89"W, 90d 0' 3.94"N)
Lower Left  (-3733200.000,-1866600.000) (  0d 0' 7.89"W, 90d 0' 3.94"S)
Upper Right ( 3733200.000, 1866600.000) (  0d 0' 7.89"E, 90d 0' 3.94"N)
Lower Right ( 3733200.000,-1866600.000) (  0d 0' 7.89"E, 90d 0' 3.94"S)
Center      (   0.0000000,   0.0000000) (180d 0' 0.00"E,  0d 0' 0.01"N)
Band 1 Block=24888x1 Type=Byte, ColorInterp=Gray
  NoData Value=0

We can create a new Custom CRS in QGIS with that WKT definition:

1

Now, when we load the raster file in QGIS, the CRS of the layer is recognized with the previously defined custom CRS (USER:100222 in my case):

2

About the map canvas CRS, I have it as 'No projection (or unknown/non-Earth projection)'. I do that because I don't want the map reprojected on-the-fly to any system, to see the coordinates of the data in the map as they are. Also, we are creating a new datum, and there are not datum transformations defined for it.


About the CSV file, we can use MASTER_LAT and MASTER_LON as geodetic coordinates for Point geometries.

We can create a new "latlon" CRS WKT, with the same geographic CRS as the mosaic, to be assigned to the CSV file to import it in QGIS:

GEOGCRS["GCS_Pluto",
    DATUM["D_Pluto",
        ELLIPSOID["Pluto",1188300,0,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001
                ]
            ]
        ]
    ],
    PRIMEM["Reference_Meridian",0,
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122
            ]
        ]
    ],
    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
        ]
    ]
]

(The WKT is based on EPSG:4326 CRS WKT, with the Geography CRS changed to match the mosaic Geography CRS.)

Now, we can create a custom geographic CRS in QGIS with that WKT:

3


We can load the CSV file and assign the new custom CRS (USER:100223 in my case) to it:

4

I am not reprojecting the canvas, and we know that the mosaic is reprojected to meters but the points are in degrees, I don't expect that they overlap. But I can see the points coordinates in the map zooming to the layer:

5


We can see that the Longitudes domain is from 0 to 360 degrees, but PROJ/GDAL/QGIS can handle them, and the center meridian of the mosaic projection is also 180 degrees. Just export the CSV to a GeoPackage with a reprojection to the Equidistant Cylindric system:

6


Now, we have a vector layer reprojected to the same Projected CRS as the mosaic, ready to work with them:

7

3
  • Hi. Many thanks for the answer! I have learned a lot with this post. But I still have a problem. Can you you zoom in to the point with MASTER_ID 690 (MASTER_LON = 130; MASTER_LAT = 7)? This image : imgur.com/rg1ieF8 is from my project and as you can see, almost all the points are still not matching the craters. I just wanto to check if I'm not missing something. – caio villaca Mar 20 '20 at 14:39
  • I see those differences. In that case, a difference of 7000m in projected X (Equidistant Cylindrical from a 1188300m radius sphere) is a difference of 0.337 degrees in longitude, but the standard deviation for that point says 0.135 degrees (although the level of confidence is 0). I also defined each CRS with the radius of the Pluto2000 datum (1195000m) and projected the points with that datum, but the differences are greater, and the mosaic has its datum defined so I did not change it. In my opinion, they are (points or mosaic) mapping errors. – Gabriel De Luca Mar 20 '20 at 15:35
  • 1
    If you can establish a relation between the differences and the coordinates, it would be possible to assume another opinion: that the error is in the datum of one or both sources. The analysis would be interesting but, at that level, I think a consultation with the data provider would be appropriate. – Gabriel De Luca Mar 20 '20 at 15:49

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.