Skip to content

Plotting in SPACE

julia
using GeoMakie, CairoMakie
using Rasters, ArchGDAL
using ImageIO, DataDeps

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

julia
DataDeps.register(
    DataDeps.DataDep(
        "vesta_image",
        """
        An image of Vesta's topography.
        Image credit: NASA/JPL-Caltech/UCLA/MPS/DLR/IDA/PSI
        """,
        "https://photojournal.jpl.nasa.gov/jpeg/PIA17037.jpg",
    )
)
DataDeps.DataDep("vesta_image", "https://photojournal.jpl.nasa.gov/jpeg/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.

julia
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).

julia
vesta_crs = GeoMakie.GeoFormatTypes.ProjString("+proj=longlat +a=285000 +b=229000 +type=crs")
GeoFormatTypes.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:

julia
ENV["PROJ_IGNORE_CELESTIAL_BODY"] = "yes"
"yes"

Now, we can create a Raster from the image.

julia
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:

julia
fig = Figure()
ga = GeoAxis(fig[1, 1]; source = vesta_crs, dest = "+proj=geos +h=35785831 +a=285000 +b=229000")
mi = meshimage!(ga, -180..180, -90..90, vesta_image_matrix; shading = NoShading)
limits!(ga, (-180, 180), (-90, 90))
ga.title = "The surface of Vesta"
fig


This page was generated using Literate.jl.