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.