Skip to content

Commit

Permalink
Improve raster example
Browse files Browse the repository at this point in the history
  • Loading branch information
asinghvi17 committed Sep 2, 2024
1 parent 1b5b936 commit e897be2
Showing 1 changed file with 35 additions and 14 deletions.
49 changes: 35 additions & 14 deletions examples/raster_warping_masking.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,12 @@ Description = "Semi involved example showing raster warping and masking, and plo
Cover = fig
```

#=
In this example, we want to show how you can create manipulate rasters for the purpose of visualization.
=#

using CairoMakie, GeoMakie
using Rasters, ArchGDAL,NaturalEarth, FileIO

Expand All @@ -25,25 +31,37 @@ blue_marble_raster = Raster(
)
# Construct a mask of land using Natural Earth land polygons
land_mask = Rasters.boolmask(NaturalEarth.naturalearth("land", 10); to = blue_marble_raster)
#
# Now, create a new raster as a copy of the background image,
# and set the values of all sea areas to be transparent.
land_raster = deepcopy(blue_marble_raster)
land_raster[(!).(land_mask)] .= RGBA{Makie.Colors.N0f8}(0.0, 0.0, 0.0, 0.0)
land_raster

# Do the same for the sea.
sea_raster = deepcopy(blue_marble_raster)
sea_raster[land_mask] .= RGBA{Makie.Colors.N0f8}(0.0, 0.0, 0.0, 0.0)
sea_raster

# ## Constructing fake data

# Now, we construct fake data.
# Now, we construct fake data.
# First, we create a matrix of values.
field = [exp(cosd(x)) + 3(y/90) for x in -180:180, y in -90:90]

# Then, we wrap it in a Raster. This is now gridded data, along the specified `X` and `Y` dimensions.
# Note the `missingval = NaN` keyword argument. This is used to indicate that if the raster does have missing values,
# they should be represented as `NaN`, and not 0, for example.
ras = Raster(field, (X(LinRange(-30, 30, size(field, 1))), Y(LinRange(-5, 5, size(field, 2)))); missingval = NaN)

# We rotate the raster and then translate it to approximate the position in the image.
rotmat = Makie.rotmatrix2d(-π/5)

# We warp the raster by the parameters below. This is in total basically
# rotating the raster by the rotation matrix `rotmat`, and then translating it
# by 30 units in the x direction and -36 in the y direction.
#
# This is where the `missingval` keyword argument comes in. If we didn't have it,
# the values that are now `NaN` because they are outside the original raster would be filled by `0`, which is not what we want.
#
# Rasters.warp is pretty much a wrapper around [GDAL's `gdalwarp` program](https://gdal.org/en/latest/programs/gdalwarp.html), so all of the options are available. It's incredibly versatile and flexible but is a bit tricky to use.
transformed = Rasters.warp(
ras,
Dict(
Expand All @@ -58,14 +76,15 @@ transformed = Rasters.warp(
)

# ## Plotting the data

fig = Figure(size = (1000, 1000))
# We use an orthographic GeoAxis to show the data.
# First, we create the figure and the axis,
fig = Figure()
ax = GeoAxis(fig[1, 1]; dest = "+proj=ortho +lon_0=0 +lat_0=0")

# Then, we plot the sea and land rasters. These are our background images.
sea_plot = meshimage!(ax, sea_raster; source = "+proj=longlat +datum=WGS84 +type=crs", npoints = 300)

land_plot = meshimage!(ax, land_raster; source = "+proj=longlat +datum=WGS84 +type=crs", npoints = 300)

# Finally, we plot the transformed data.
data_plot = surface!(ax, transformed; shading = NoShading)
# Get the land plot above all other plots in the geoaxis.
translate!(land_plot, 0, 0, 1)
Expand All @@ -74,13 +93,15 @@ fig

# For extra credit, we'll find the center of the raster and transform the orthographic projection to be centered on it:

using Statistics
xy_matrix = Rasters.DimensionalData.DimPoints(transformed) |> collect
center_x = mean(first.(xy_matrix))
center_y = mean(last.(xy_matrix))

# We'll use the `ortho` projection, which is centered on the point we just found.
ax.dest = "+proj=ortho +lon_0=$(center_x) +lat_0=$(center_y)"

ax.dest[] = "+proj=ortho +lon_0=$(center_x) +lat_0=$(center_y)"
# Somehow, the limits are set to the wrong thing, so we'll set them manually and let the axis take care of it.
ax.limits[] = (-180, 180, -90, 90)
fig

# ## Plotting on the 3D globe
Expand All @@ -90,20 +111,20 @@ fig
using Geodesy

# fig = Figure(size = (1000, 1000));
fig = Figure()
fig2 = Figure()

ax = LScene(fig[1, 1])
ax = LScene(fig2[1, 1])

sea_plot = meshimage!(ax, sea_raster; npoints = 300)

land_plot = meshimage!(ax, land_raster; npoints = 300)

data_plot = surface!(ax, transformed; shading = NoShading)

land_plot.z_level[] = 1000
land_plot.z_level[] = 1000 # this means raising land by 1 km over the sea!

sea_plot.transformation.transform_func[] = Geodesy.ECEFfromLLA(Geodesy.WGS84())
land_plot.transformation.transform_func[] = Geodesy.ECEFfromLLA(Geodesy.WGS84())
data_plot.transformation.transform_func[] = Geodesy.ECEFfromLLA(Geodesy.WGS84())

fig
fig2

0 comments on commit e897be2

Please sign in to comment.