Plotting in SPACE!

This example shows how to plot a 2D image of Vesta's surface in a 3D globe. It uses the GlobeAxis type to create a globe axis and plot the image on it. The ellipsoid parameters are set to the correct values for Vesta, which is a pretty oblate ellipsoid.

using GeoMakie, GLMakie
using Rasters
using ImageIO, DataDeps

We use the DataDeps.jl package to register the image, and avoid repeated downloads - especially on CI.

DataDeps.register(
    DataDeps.DataDep(
        "vesta_image",
        """
        An image of Vesta's topography.
        Image credit: NASA/JPL-Caltech/UCLA/MPS/DLR/IDA/PSI
        """,
        "https://assets.science.nasa.gov/content/dam/science/psd/photojournal/pia/pia17/pia17037/PIA17037.jpg",
    )
)

DataDeps.DataDep("vesta_image", "https://assets.science.nasa.gov/content/dam/science/psd/photojournal/pia/pia17/pia17037/PIA17037.jpg", nothing, DataDeps.fetch_default, identity, "An image of Vesta's topography.\nImage credit: NASA/JPL-Caltech/UCLA/MPS/DLR/IDA/PSI\n")

Now, we can load the image and create a Raster from it. This is stored in memory.

vesta_image_matrix = FileIO.load(joinpath(DataDeps.datadep"vesta_image", "PIA17037.jpg"))

Since this is an image of Vesta, we need to define its coordinate reference system (CRS). We use Proj's ellipsoid parameters to define the ellipsoid of Vesta, specifically the semi-major and semi-minor axes (+a and +b respectively).

vesta_crs = GeoMakie.GeoFormatTypes.ProjString("+proj=longlat +a=285000 +b=229000 +type=crs")

ProjString: +proj=longlat +a=285000 +b=229000 +type=crs

Additionally, since the Vesta CRS is on a non-Earth datum, we have to set this environment variable so that Proj knows that we are aware of this problem:

ENV["PROJ_IGNORE_CELESTIAL_BODY"] = "yes"

"yes"

Now, we can create a Raster from the image.

vesta_raster = Raster(
    rotr90(vesta_image_matrix);
    dims = (
        X(LinRange(-180, 180, size(vesta_image_matrix, 2))),
        Y(LinRange(-90, 90, size(vesta_image_matrix, 1)))
    ),
    crs = vesta_crs
)

Now that we have the data, we can go ahead and plot it in a GlobeAxis with the right ellipsoid params!

f, a, p = meshimage(
    -180..180, -90..90, vesta_image_matrix;
    uv_transform = :rotr90,
    axis = (;
        type = GlobeAxis,
        dest = GeoMakie.Geodesy.Ellipsoid(; a = "285000", b="229000"),
        camera_longlat = (0, 0),
    )
)

---

This page was generated using Literate.jl.