Skip to content

Commit

Permalink
resample examples, fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@lazarusA
lazarusA committed Dec 30, 2024
1 parent f441146 commit 68410c4
Show file tree
Hide file tree
Showing 2 changed files with 188 additions and 53 deletions.
57 changes: 51 additions & 6 deletions docs/src/proje.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,54 @@ using DimensionalData.Lookups
using NaNStatistics
using GLMakie

# resample section
ras = Raster(WorldClim{BioClim}, 5)
ras_m = replace_missing(ras, missingval=NaN)
SINUSOIDAL_CRS = ProjString("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
ras_m = replace_missing(ras, missingval=NaN) # ArchGDAL warp handles NaNs via the metadata in the Raster

ras_sample = resample(ras_m; size=(2160, 1080))

methods = ["average", "mode", "max", "sum"]
sizes = [(2160, 1080), (1440, 720), (720, 360), (360, 180)]
resolutions = [0.16666666666666666, 0.25, 0.5, 1.0]

# other available methods to try: "bilinear", "cubic", "cubicspline", "lanczos", "min", "med", "q1", "q3", "near".

method_sizes = [resample(ras_m; size=size, method=method) for method in methods for size in sizes]
method_res = [resample(ras_m; res=res, method=method) for method in methods for res in resolutions]

with_theme(Rasters.theme_rasters()) do
colorrange = (nanminimum(ras_m), nanmaximum(ras_m))
hm=nothing
fig = Figure(; size = (1000, 600))
axs = [Axis(fig[i,j], title="size=$(size), method=:$(method)", titlefont=:regular)
for (i, method) in enumerate(methods) for (j, size) in enumerate(sizes)]
for (i, ax) in enumerate(axs)
hm = heatmap!(ax, method_sizes[i]; colorrange)
end
Colorbar(fig[:,end+1], hm)
hidedecorations!.(axs; grid=false)
rowgap!(fig.layout, 5)
colgap!(fig.layout, 10)
fig
end

with_theme(Rasters.theme_rasters()) do
colorrange = (nanminimum(ras_m), nanmaximum(ras_m))
hm=nothing
fig = Figure(; size = (1000, 600))
axs = [Axis(fig[i,j], title="res=$(round(res, digits=4)), method=:$(method)", titlefont=:regular)
for (i, method) in enumerate(methods) for (j, res) in enumerate(resolutions)]
for (i, ax) in enumerate(axs)
hm = heatmap!(ax, method_res[i]; colorrange)
end
Colorbar(fig[:,end+1], hm)
hidedecorations!.(axs; grid=false)
rowgap!(fig.layout, 5)
colgap!(fig.layout, 10)
fig
end

# SINUSOIDAL_CRS = ProjString("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")

SINUSOIDAL_CRS = ProjString("+proj=sinu +lon_0=0 +type=crs")

Expand All @@ -19,12 +64,13 @@ ras_sin = resample(ras_m; size=(2160, 1080), crs=SINUSOIDAL_CRS, method="average

fig, ax, plt = heatmap(ras_sin)
Colorbar(fig[1,2], plt)
fig
display(GLMakie.Screen(), fig)


ras_epsg = resample(ras_sin; size=(1440,720), crs=EPSG(4326), method="average")

heatmap(ras_epsg)
fig = heatmap(ras_epsg)
display(GLMakie.Screen(), fig)

locus_resampled = DimensionalData.shiftlocus(Center(), ras_epsg)
heatmap(locus_resampled)
Expand All @@ -34,8 +80,7 @@ y_range = LinRange(89.75, -90, 720)
ras_data = ras_epsg.data

ras_scratch = Raster(ras_data, (X(x_range; sampling=Intervals(Start())),
Y(y_range; sampling=Intervals(Start()))), crs=EPSG(4326))

Y(y_range; sampling=Intervals(Start()))), crs=EPSG(4326), missingval=NaN)
heatmap(ras_scratch)
# locus_resampled = DimensionalData.shiftlocus(Center(), ras_scratch)

Expand Down
184 changes: 137 additions & 47 deletions docs/src/tutorials/resample_warp.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,19 +3,7 @@ CollapsedDocStrings=true
```
# Reprojection and resampling

- What is resampling?
- When to resample vs reproject
- Things to keep in mind
- GDAL always changes the locus to cell sampling, you can reset this by using `shiftlocus`
- You can in fact resample to another raster, if you want perfect alignment.
- This doesn't work for irregularly sampled rasters.
- Resampling is a lossy operation and takes time. Try to avoid repeatedly resampling, and if you must do so, crop or trim the raster as much as you can first.
- Show the different resampling methods, in a grid
- Show some different projections and ways of constructing them
- Show how to use `size` and `res` to change the resolution of a raster
- Show how to use `warp` to reproject a raster

## What is resampling?
### What is resampling?

**[`resample`](@ref)** "re-samples" the
data by interpolation and can also aggregate or disaggregate, changing the resolution.
Expand All @@ -25,6 +13,12 @@ resampling methods.
This uses GDAL's `gdalwarp` algorithm under the hood. You can call that via [`warp`](@ref)
if you need more control, but generally `resample` is sufficient.

::: tip warp, contributions are welcome!

- Show how to use `warp` to reproject a raster

:::

Rasters.jl has a few other methods to change the lookups of a raster. These are:
- [`reproject`](@ref), which directly reprojects the lookup axes
(but is **only usable for specific cases**, where the source and destination
Expand All @@ -40,7 +34,7 @@ Rasters.jl has a few other methods to change the lookups of a raster. These are
Of all these methods, **`resample`** is the most flexible and powerful, and is what we will focus on here.
It is, however, also the slowest. So if another method is applicable to your problem, you should consider it.

## How `resample` works
### How `resample` works

`resample` uses GDAL's `gdalwarp` algorithm under the hood. This is a battle-tested algorithm
and is generally pretty robust. However, it has the following limitations:
Expand All @@ -49,95 +43,191 @@ and is generally pretty robust. However, it has the following limitations:
- It can only accept some primitive types for the input data, since that data is passed directly to a C library.
Things like `RGB` or user-defined types are not usually supported.

`resample` allows you to specify several methods, which you can see if you expand the docstring below.
`resample` allows you to specify several methods, see some of them in the next section.

## `resample` and projections with `ProjString`
### `resolution`, `size` and `methods`

Geospatial datasets will come in different [projections](https://proj.org/en/9.4/operations/projections/index.html) or coordinate reference systems (CRS) for many reasons. Here, we will focus on `MODIS SINUSOIDAL` and `EPSG`, and transformations between them.

Let's start by loading the neccesary packages:
Let's start by loading the necessary packages:

````@example modis
````@example resample
using Rasters, RasterDataSources, ArchGDAL
using DimensionalData
using DimensionalData.Lookups
using NaNStatistics
using CairoMakie
````

and let's load a test raster

````@example modis
````@example resample
ras = Raster(WorldClim{BioClim}, 5)
ras_m = replace_missing(ras, missingval=NaN)
````

and let's also take a look
resampling to a given `size` or `res ≡ resolution` providing a `method` is done with

````@example modis
fig, ax, plt = heatmap(ras_m)
Colorbar(fig[1,2], plt)
fig
:::tabs

== tab size

````@ansi resample
ras_sample = resample(ras_m; size=(2160, 1080), method="average")
````

### MODIS SINUSOIDAL PROJECTION
== tab resolution

````@example modis
# ? is this the right ProjString ?
# SINUSOIDAL_CRS = ProjString("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
````@ansi resample
ras_sample = resample(ras_m; res=1.0, method="average")
````

:::

other available methods to try: `"mode"`, `"max"`, `"sum"`, `"bilinear"`, `"cubic"`, `"cubicspline"`, `"lanczos"`, `"min"`, `"med"`, `"q1"`, `"q3"` and `"near"`.

Let's consider a few more examples, with the following options:

````@example resample
methods = ["average", "mode", "max", "sum"]
sizes = [(2160, 1080), (1440, 720), (720, 360), (360, 180)]
resolutions = [0.16666666666666666, 0.25, 0.5, 1.0];
nothing # hide
````

:::tabs

== tab sizes and methods

````@example resample
method_sizes = [resample(ras_m; size=size, method=method) for method in methods for size in sizes]
with_theme(Rasters.theme_rasters()) do
colorrange = (nanminimum(ras_m), nanmaximum(ras_m))
hm=nothing
fig = Figure(; size = (1000, 600))
axs = [Axis(fig[i,j], title="size=$(size), method=:$(method)", titlefont=:regular)
for (i, method) in enumerate(methods) for (j, size) in enumerate(sizes)]
for (i, ax) in enumerate(axs)
hm = heatmap!(ax, method_sizes[i]; colorrange)
end
Colorbar(fig[:,end+1], hm)
hidedecorations!.(axs; grid=false)
rowgap!(fig.layout, 5)
colgap!(fig.layout, 10)
fig
end
````

== tab resolutions and methods

````@example resample
method_res = [resample(ras_m; res=res, method=method) for method in methods for res in resolutions]
with_theme(Rasters.theme_rasters()) do
colorrange = (nanminimum(ras_m), nanmaximum(ras_m))
hm=nothing
fig = Figure(; size = (1000, 600))
axs = [Axis(fig[i,j], title="res=$(round(res, digits=4)), method=:$(method)", titlefont=:regular)
for (i, method) in enumerate(methods) for (j, res) in enumerate(resolutions)]
for (i, ax) in enumerate(axs)
hm = heatmap!(ax, method_res[i]; colorrange)
end
Colorbar(fig[:,end+1], hm)
hidedecorations!.(axs; grid=false)
rowgap!(fig.layout, 5)
colgap!(fig.layout, 10)
fig
end
````

:::


## `reproject` with `resample` using a `ProjString`

Geospatial datasets will come in different [projections](https://proj.org/en/9.4/operations/projections/index.html) or coordinate reference systems (CRS) for many reasons. Here, we will focus on `MODIS SINUSOIDAL` and `EPSG`, and transformations between them.

Let's load our test raster

````@example resample
ras = Raster(WorldClim{BioClim}, 5)
ras_m = replace_missing(ras, missingval=NaN);
nothing # hide
````

### Sinusoidal Projection (MODIS)

````@example resample
SINUSOIDAL_CRS = ProjString("+proj=sinu +lon_0=0 +type=crs")
````

and hence the `resample` is performed with
::: details Raw MODIS ProjString

````julia
SINUSOIDAL_CRS = ProjString("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
````

:::


and the `resample` is performed with

````@example modis
````@example resample
ras_sin = resample(ras_m; size=(2160, 1080), crs=SINUSOIDAL_CRS, method="average") # ? do again, failing locally on first try
ras_sin = resample(ras_m; size=(2160, 1080), crs=SINUSOIDAL_CRS, method="average") # hide
nothing # hide
````

::: tip

`GDAL` always changes the locus to cell sampling, you can reset this by using `shiftlocus`.

:::

let's compare the total counts!

````@example modis
````@example resample
nansum(ras_m), nansum(ras_sin)
````

and, how does this looks like?

````@example modis
````@example resample
fig, ax, plt = heatmap(ras_sin)
Colorbar(fig[1,2], plt)
fig
````

now, let's go back to `latitude` and `longitude` and reduce the resolution

````@example modis
````@example resample
ras_epsg = resample(ras_sin; size=(1440,720), crs=EPSG(4326), method="average")
````

and let's apply `shiftlocus` such that the lookups share the exact same grid, which might be needed when building bigger datasets:

````@example modis
````@example resample
locus_resampled = DimensionalData.shiftlocus(Center(), ras_epsg)
````

and compare the total counts!

````@example modis
````@example resample
nansum(ras_m), nansum(locus_resampled)
````

````@example modis
::: info Things to keep in mind

- You can in fact resample to another raster `resample(ras; to=ref_ras)`, if you want perfect alignment. Contributions are welcome for this use case!
- This doesn't work for irregularly sampled rasters.

:::


````@example resample
fig, ax, plt = heatmap(ras_epsg)
Colorbar(fig[1,2], plt)
fig
````

### Construct a Raster from scratch natively in the sinusoidal projection
### A `Raster` from scratch

````@example modis
````@example resample
x_range = LinRange(-180, 179.75, 1440)
y_range = LinRange(89.75, -90, 720)
ras_data = ras_epsg.data
Expand All @@ -148,7 +238,7 @@ create the raster

````@ansi modis
ras_scratch = Raster(ras_data, (X(x_range; sampling=Intervals(Start())),
Y(y_range; sampling=Intervals(Start()))), crs=EPSG(4326))
Y(y_range; sampling=Intervals(Start()))), crs=EPSG(4326), missingval=NaN)
````

Expand All @@ -162,7 +252,7 @@ This requires that you run `using Rasters.Lookups`, where the `Intervals` and `S

and take a look

````@example modis
````@example resample
fig, ax, plt = heatmap(ras_scratch)
Colorbar(fig[1,2], plt)
fig
Expand All @@ -174,15 +264,15 @@ and the corresponding resampled projection
ras_sin_s = resample(ras_scratch; size=(1440,720), crs=SINUSOIDAL_CRS, method="average")
````

````@example modis
````@example resample
fig, ax, plt = heatmap(ras_sin_s)
Colorbar(fig[1,2], plt)
fig
````

and go back from `sin` to `epsg`:

````@example modis
````@example resample
ras_epsg = resample(ras_sin_s; size=(1440,720), crs=EPSG(4326), method="average")
locus_resampled = DimensionalData.shiftlocus(Center(), ras_epsg)
Expand All @@ -194,7 +284,7 @@ fig

and compare the total counts again!

````@example modis
````@example resample
nansum(ras_m), nansum(locus_resampled)
````

Expand Down

0 comments on commit 68410c4

Please sign in to comment.