From 6e91cb9e06420bb8c2129caadd1e0882288c2be5 Mon Sep 17 00:00:00 2001 From: Lazaro Alonso Date: Sun, 29 Dec 2024 17:03:50 +0100 Subject: [PATCH] more orga --- docs/src/.vitepress/config.mts | 6 + docs/src/manual/cellarea.md | 38 +++--- docs/src/manual/data_sources.md | 4 - docs/src/manual/resample_warp.md | 67 ---------- docs/src/resample_proj.md | 126 ------------------- docs/src/tutorials/resample_warp.md | 184 ++++++++++++++++++++++++++++ src/extensions.jl | 11 +- 7 files changed, 215 insertions(+), 221 deletions(-) delete mode 100644 docs/src/manual/resample_warp.md delete mode 100644 docs/src/resample_proj.md create mode 100644 docs/src/tutorials/resample_warp.md diff --git a/docs/src/.vitepress/config.mts b/docs/src/.vitepress/config.mts index f26f08ff1..b0eaeb13d 100644 --- a/docs/src/.vitepress/config.mts +++ b/docs/src/.vitepress/config.mts @@ -19,6 +19,7 @@ const navTemp = { { text: 'Manual', items: [ { text: 'Methods', link: '/manual/methods' }, + { text: 'cellarea', link: '/manual/cellarea' }, { text: 'Array Operations', link: '/manual/array_operations' }, { text: 'Data Sources', link: '/manual/data_sources' }, { text: 'Plots', @@ -31,6 +32,8 @@ const navTemp = { }, { text: 'Tutorials', items: [ + { text: 'Spatial mean', link: '/tutorials/spatial_mean' }, + { text: 'Reprojection and resampling', link: '/tutorials/resample_warp'}, { text: 'Species Distribution Modelling', link: '/tutorials/gbif_wflow' }, ] }, @@ -99,6 +102,7 @@ export default defineConfig({ { text: 'Manual', items: [ { text: 'Methods', link: '/manual/methods' }, + { text: 'cellarea', link: '/manual/cellarea' }, { text: 'Array Operations', link: '/manual/array_operations' }, { text: 'Data Sources', link: '/manual/data_sources' }, { text: 'Plots', @@ -111,6 +115,8 @@ export default defineConfig({ }, { text: 'Tutorials', items: [ + { text: 'Spatial mean', link: '/tutorials/spatial_mean' }, + { text: 'Reprojection and resampling', link: '/tutorials/resample_warp'}, { text: 'Species Distribution Modelling', link: '/tutorials/gbif_wflow' }, ] }, diff --git a/docs/src/manual/cellarea.md b/docs/src/manual/cellarea.md index 6a14df3bb..a1a2475d0 100644 --- a/docs/src/manual/cellarea.md +++ b/docs/src/manual/cellarea.md @@ -8,29 +8,24 @@ CollapsedDocStrings=true cellarea ``` -`cellarea` computes the area of each cell in a raster. -This is useful for a number of reasons - if you have a variable like +Computing the area of each cell in a raster is useful for a number of reasons - if you have a variable like population per cell, or elevation ([spatially extensive variables](https://r-spatial.org/book/05-Attributes.html#sec-extensiveintensive)), you'll want to account for the fact that different cells have different areas. -`cellarea` returns a Raster with the same x and y dimensions as the input, -where each value in the raster encodes the area of the cell (in meters by default). - -You can specify whether you want to compute the area in the plane of your projection -(`Planar()`), on a sphere of some radius (`Spherical(; radius=...)`). - - Let's construct a raster and see what this looks like! We'll keep it in memory. -The spherical method relies on the [Proj.jl](https://github.com/JuliaGeo/Proj.jl) package to perform coordinate transformation, so that has to be loaded explicitly. +The spherical method relies on the [Proj.jl](https://github.com/JuliaGeo/Proj.jl) package to perform coordinate transformation, so that has to be loaded explicitly. -````@example cella +````@example cellarea using Rasters, Proj ```` -To construct a raster, we'll need to specify the x and y dimensions. These are called "lookups" in Rasters.jl. +To construct a raster, we'll need to specify the `x` and `y` dimensions. These are called `lookups` in `Rasters.jl.` + +````@example cellarea using Rasters.Lookups +```` + We can now construct the x and y lookups. Here we'll use a start-at-one, step-by-five grid. Note that we're specifying that the "sampling", i.e., what the coordinates actually mean, is `Intervals(Start())`, meaning that the coordinates are the starting point of each interval. @@ -38,33 +33,34 @@ is `Intervals(Start())`, meaning that the coordinates are the starting point of This is in contrast to `Points()` sampling, where each index in the raster represents the value at a sampling point; here, each index represents a grid cell, which is defined by the coordinate being at the start. -````@example cella +````@example cellarea x = X(1:5:30; sampling = Intervals(Start()), crs = EPSG(4326)) -y = Y(50:5:80; sampling = Intervals(Start()), crs = EPSG(4326)) +y = Y(50:5:80; sampling = Intervals(Start()), crs = EPSG(4326)); +nothing # hide ```` I have chosen the y-range here specifically so we can show the difference between spherical and planar `cellarea`. -````@example cella +````@ansi cellarea ras = Raster(ones(x, y); crs = EPSG(4326)) ```` We can just call `cellarea` on this raster, which returns cell areas in meters, on Earth, assuming it's a sphere: -````@example cella +````@ansi cellarea cellarea(ras) ```` and if we plot it, you can see the difference in cell area as we go from the equator to the poles: -````@example cella -using Makie, CairoMakie# , GeoMakie -heatmap(ras; axis = (; aspect = DataAspect())) +````@example cellarea +using CairoMakie +heatmap(cellarea(ras); axis = (; aspect = DataAspect())) ```` We can also try this using the planar method, which simply computes the area of the rectangle using `area = x_side_length * y_side_length`: -````@example cella +````@ansi cellarea cellarea(Planar(), ras) ```` diff --git a/docs/src/manual/data_sources.md b/docs/src/manual/data_sources.md index 2604d9ca0..93df91881 100644 --- a/docs/src/manual/data_sources.md +++ b/docs/src/manual/data_sources.md @@ -43,10 +43,6 @@ They are always 3 dimensional, and have `Y`, `X` and [`Band`](@ref) dimensions. using Rasters ```` -````@docs -smapseries -```` - ## Writing file formats to disk Files can be written to disk with ArchGDAL.jl and NCDatasets.jl backends using diff --git a/docs/src/manual/resample_warp.md b/docs/src/manual/resample_warp.md deleted file mode 100644 index 728af85b6..000000000 --- a/docs/src/manual/resample_warp.md +++ /dev/null @@ -1,67 +0,0 @@ -## `resample` warping a raster - -```@meta -CollapsedDocStrings=true -``` - -````@example resample -using Rasters, ArchGDAL -using RasterDataSources -using NaNStatistics -using CairoMakie -```` - -## What is resampling? - -**[`resample`](@ref)** "re-samples" the -data by interpolation and can also aggregate or disaggregate, changing the resolution. -It always returns a `Regular` lookup (like a range), and is the most flexible of the -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. - -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 - coordinate systems are both cylindrical, like the long-lat, Mercator, or Web-Mercator projections.) - - This is a lossless operation and keeps the data exactly the same - only the axes are changed. - -- [`aggregate`](@ref) and [`disaggregate`](@ref), which change the resolution of - the raster by merging ([`aggregate`](@ref)) or splitting ([`disaggregate`](@ref)) cells. - - They can't change cells fractionally, and can't change the projection or coordinate system. - -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 - -`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: -- It always assumes cell-based sampling, instead of point-based sampling. This does mean that - point-based rasters are converted to cell-based sampling. -- 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. - -```@docs; canonical=false -resample -``` - -Topics: -- 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 - - diff --git a/docs/src/resample_proj.md b/docs/src/resample_proj.md deleted file mode 100644 index 544e63ddf..000000000 --- a/docs/src/resample_proj.md +++ /dev/null @@ -1,126 +0,0 @@ -## `resample` and projections with `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 start by loading the neccesary packages: - -````@example modis -using Rasters, RasterDataSources, ArchGDAL -using DimensionalData -using DimensionalData.Lookups -using NaNStatistics -using CairoMakie -```` - -and let's load a test raster - -````@example modis -ras = Raster(WorldClim{BioClim}, 5) -ras_m = replace_missing(ras, missingval=NaN) -```` - -and let's also take a look - -````@example modis -fig = plot(ras_m; colorrange=(0,100)) -```` - -### MODIS SINUSOIDAL PROJECTION - -````@example -# ? is this the right ProjString ?, do we need to shift lat, lon? -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 hence the `resample` is performed with - -````@example modis -ras_sin = resample(ras_m; size=(1440,720), crs=SINUSOIDAL_CRS, method="sum") -```` - -let's compare the total counts! - -````@example modis -nansum(ras_m), nansum(ras_sin) -```` - -and, how does this looks like? - -````@example modis -fig = plot(ras_sin; colorrange=(0,100)) -```` - -now, let's go back to `latitude` and `longitude` - -````@example modis -# ? also here, do we need to shift X, Y? -ras_epsg = resample(ras_sin; size=(1440,720), crs=EPSG(4326), method="sum") -```` - -and let's apply `shiftlocus` such that the lookups share the exact same grid, which might be needed when building bigger datasets: - -````@example modis -locus_resampled = DimensionalData.shiftlocus(Center(), ras_epsg) -```` - -and compare the total counts! - -````@example modis -nansum(ras_m), nansum(locus_resampled) -```` - -````@example modis -fig = plot(ras_epsg; colorrange=(0,100)) -```` - -### Construct a Raster from scratch natively in the sinusoidal projection - -````@example modis -x_range = range(-2.0015109355797417e7, 1.998725401355172e7, 1440) -y_range = range(9.979756529777847e6, -1.0007554677898709, 720) -ra_data = ras_sin.data; -nothing # hide -```` - -create the raster - -````@example modis -ras_scratch = Raster(ra_data, (X(x_range; sampling=Intervals(Start())), Y(y_range; sampling=Intervals(Start()))), - crs=SINUSOIDAL_CRS) -```` - -::: warning -At the moment, you need to specify `sampling=Intervals(Start())` for `X` and `Y`. - -This requires that you run `using Rasters.Lookups`, where the `Intervals` and `Start` types are defined. -::: - -and take a look - -````@example modis -fig = plot(ras_scratch; colorrange=(0,100)) -```` - -and the corresponding resampled projection - -````@example modis -ras_latlon = resample(ras_scratch; size=(1440,720), crs=EPSG(4326), method="sum") -locus_resampled = DimensionalData.shiftlocus(Center(), ras_latlon) -```` - -````@example modis -fig = plot(ras_latlon; colorrange=(0,100)) -```` - -and compare the total counts again! - -````@example modis -nansum(ras_m), nansum(locus_resampled) -```` - -::: danger -Note that all counts are a little bit off. Could we mitigate this some more? -::: - - - diff --git a/docs/src/tutorials/resample_warp.md b/docs/src/tutorials/resample_warp.md new file mode 100644 index 000000000..c4e8c4d2f --- /dev/null +++ b/docs/src/tutorials/resample_warp.md @@ -0,0 +1,184 @@ +```@meta +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? + +**[`resample`](@ref)** "re-samples" the +data by interpolation and can also aggregate or disaggregate, changing the resolution. +It always returns a `Regular` lookup (like a range), and is the most flexible of the +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. + +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 + coordinate systems are both cylindrical, like the long-lat, Mercator, or Web-Mercator projections.) + + This is a lossless operation and keeps the data exactly the same - only the axes are changed. + +- [`aggregate`](@ref) and [`disaggregate`](@ref), which change the resolution of + the raster by merging ([`aggregate`](@ref)) or splitting ([`disaggregate`](@ref)) cells. + + They can't change cells fractionally, and can't change the projection or coordinate system. + +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 + +`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: +- It always assumes cell-based sampling, instead of point-based sampling. This does mean that + point-based rasters are converted to cell-based sampling. +- 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` and projections with `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 start by loading the neccesary packages: + +````@example modis +using Rasters, RasterDataSources, ArchGDAL +using DimensionalData +using DimensionalData.Lookups +using NaNStatistics +using CairoMakie +```` + +and let's load a test raster + +````@example modis +ras = Raster(WorldClim{BioClim}, 5) +ras_m = replace_missing(ras, missingval=NaN) +```` + +and let's also take a look + +````@example modis +fig = plot(ras_m; colorrange=(0,100)) +```` + +### MODIS SINUSOIDAL PROJECTION + +````@example modis +# ? is this the right ProjString ?, do we need to shift lat, lon? +# 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") +```` + +and hence the `resample` is performed with + +````@example modis +ras_sin = resample(ras_m; size=(1440,720), crs=SINUSOIDAL_CRS, method="sum") +```` + +let's compare the total counts! + +````@example modis +nansum(ras_m), nansum(ras_sin) +```` + +and, how does this looks like? + +````@example modis +fig = plot(ras_sin; colorrange=(0,100)) +```` + +now, let's go back to `latitude` and `longitude` + +````@example modis +# ? also here, do we need to shift X, Y? +ras_epsg = resample(ras_sin; size=(1440,720), crs=EPSG(4326), method="sum") +```` + +and let's apply `shiftlocus` such that the lookups share the exact same grid, which might be needed when building bigger datasets: + +````@example modis +locus_resampled = DimensionalData.shiftlocus(Center(), ras_epsg) +```` + +and compare the total counts! + +````@example modis +nansum(ras_m), nansum(locus_resampled) +```` + +````@example modis +fig = plot(ras_epsg; colorrange=(0,100)) +```` + +### Construct a Raster from scratch natively in the sinusoidal projection + +````@example modis +x_range = range(-2.0015109355797417e7, 1.998725401355172e7, 1440) +y_range = range(9.979756529777847e6, -1.0007554677898709, 720) +ra_data = ras_sin.data; +nothing # hide +```` + +create the raster + +````@example modis +ras_scratch = Raster(ra_data, (X(x_range; sampling=Intervals(Start())), Y(y_range; sampling=Intervals(Start()))), + crs=SINUSOIDAL_CRS) +```` + +::: warning + +At the moment, you need to specify `sampling=Intervals(Start())` for `X` and `Y`. + +This requires that you run `using Rasters.Lookups`, where the `Intervals` and `Start` types are defined. + +::: + +and take a look + +````@example modis +fig = plot(ras_scratch; colorrange=(0,100)) +```` + +and the corresponding resampled projection + +````@example modis +ras_latlon = resample(ras_scratch; size=(1440,720), crs=EPSG(4326), method="sum") +locus_resampled = DimensionalData.shiftlocus(Center(), ras_latlon) +```` + +````@example modis +fig = plot(ras_latlon; colorrange=(0,100)) +```` + +and compare the total counts again! + +````@example modis +nansum(ras_m), nansum(locus_resampled) +```` + +::: danger + +Note that all counts are a little bit off. Could we mitigate this some more? + +::: + + + diff --git a/src/extensions.jl b/src/extensions.jl index f90a5ba9f..30f31fdd1 100644 --- a/src/extensions.jl +++ b/src/extensions.jl @@ -162,14 +162,19 @@ warp(args...; kw...) = throw_extension_error(warp, "ArchGDAL", :RastersArchGDALE cellarea([method], x) Gives the approximate area of each gridcell of `x`. -By assuming the earth is a sphere, it approximates the true size to about 0.1%, depending on latitude. +By assuming the earth is a sphere, it approximates the true size to about 0.1%, depending on latitude. -Run `using ArchGDAL` to make this method fully available. +Run `using ArchGDAL` or `using Proj` to make this method fully available. + +- `method`: You can specify whether you want to compute the area in the plane of your projection `Planar()` or on a sphere of some radius `Spherical(; radius=...)`(the default). -`method` can be `Spherical(; radius)` (the default) or `Planar()`. - `Spherical` will compute cell area on the sphere, by transforming all points back to long-lat. You can specify the radius by the `radius` keyword argument here. By default, this is `6371008.8`, the mean radius of the Earth. + - `Planar` will compute cell area in the plane of the CRS you have chosen. Be warned that this will likely be incorrect for non-equal-area projections. +Returns a Raster with the same x and y dimensions as the input, +where each value in the raster encodes the area of the cell (in meters by default). + ## Example ```julia