Skip to content

Commit

Permalink
dont do src files
Browse files Browse the repository at this point in the history
  • Loading branch information
lazarusA committed Dec 28, 2024
1 parent d606851 commit 8ca0e24
Show file tree
Hide file tree
Showing 2 changed files with 10 additions and 66 deletions.
64 changes: 5 additions & 59 deletions src/extensions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -164,21 +164,18 @@ warp(args...; kw...) = throw_extension_error(warp, "ArchGDAL", :RastersArchGDALE
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.
Run `using Proj` to make this method fully available.
Run `using ArchGDAL` to make this method fully available.
`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.
## Example
```jldoctest cellarea
using Rasters, Proj, Rasters.Lookups
# First, we construct a raster with cell-based sampling (`Intervals`),
# where the value represents the start of the cell in each dimension (`Start`).
xdim = X(Projected(90.0:10.0:120; sampling=Intervals(Start()), crs=EPSG(4326))) # construct cell-based sampling,
```julia
using Rasters, ArchGDAL, Rasters.Lookups
xdim = X(Projected(90.0:10.0:120; sampling=Intervals(Start()), crs=EPSG(4326)))
ydim = Y(Projected(0.0:10.0:50; sampling=Intervals(Start()), crs=EPSG(4326)))
# Now, we construct a raster from these dimensions, and compute the cell area.
myraster = rand(xdim, ydim)
cs = cellarea(myraster)
Expand All @@ -199,57 +196,6 @@ cs = cellarea(myraster)
110.0 1.23017e6 1.19279e6 1.11917e6 1.01154e6 873182.0 708290.0
120.0 1.23017e6 1.19279e6 1.11917e6 1.01154e6 873182.0 708290.0
```
Note how the cell area decreases with increasing latitude, i.e., as we move away from the equator.
Compare this to the output from `cellarea(Planar(), ...)`:
```jldoctest cellarea
cs = cellarea(Planar(), myraster)
# output
╭───────────────────────╮
│ 4×6 Raster{Float64,2} │
├───────────────────────┴────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── dims ┐
↓ X Projected{Float64} 90.0:10.0:120.0 ForwardOrdered Regular Intervals{Start},
→ Y Projected{Float64} 0.0:10.0:50.0 ForwardOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── raster ┤
extent: Extent(X = (90.0, 130.0), Y = (0.0, 60.0))
crs: EPSG:4326
└────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
↓ → 0.0 10.0 20.0 30.0 40.0 50.0
90.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0
110.0 100.0 100.0 100.0 100.0 100.0 100.0
120.0 100.0 100.0 100.0 100.0 100.0 100.0
```
The output here is clearly incorrect - but you can see how it's calculated, and that in an equal-area projection, it could be useful.
Let's look at `cellarea` in the Web Mercator projection:
```jldoctest cellarea
mywebmerc = reproject(myraster; crs = EPSG(3857)) # this is the EPSG code for Web Mercator
cs = cellarea(mywebmerc)
# output
╭───────────────────────╮
│ 4×6 Raster{Float64,2} │
├───────────────────────┴────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── dims ┐
↓ X Projected{Float64} [1.0018754171394622e7, 1.1131949079327356e7, 1.2245143987260092e7, 1.3358338895192828e7] ForwardOrdered Regular Intervals{Start},
→ Y Projected{Float64} [0.0, 1.1188899748579594e6, …, 4.865942279503176e6, 6.446275841017159e6] ForwardOrdered Irregular Intervals{Start}
├──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── raster ┤
extent: Extent(X = (1.0018754171394622e7, 1.4471533803125564e7), Y = (0.0, 8.399737889818357e6))
crs: EPSG:3857
└────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
↓ → 0.0 1.11889e6 2.27303e6 3.50355e6 4.86594e6 6.44628e6
1.00188e7 1.2332e12 1.1952e12 1.12048e12 1.01158e12 8.72085e11 7.06488e11
1.11319e7 1.2332e12 1.1952e12 1.12048e12 1.01158e12 8.72085e11 7.06488e11
1.22451e7 1.2332e12 1.1952e12 1.12048e12 1.01158e12 8.72085e11 7.06488e11
1.33583e7 1.2332e12 1.1952e12 1.12048e12 1.01158e12 8.72085e11 7.06488e11
```
$EXPERIMENTAL
"""
cellarea(x; kw...) = cellarea(Spherical(), x; kw...)
Expand Down Expand Up @@ -352,4 +298,4 @@ struct AffineProjected{T,F,A<:AbstractVector{T},M,C,MC,P,D} <: LA.Unaligned{T,1}
mappedcrs::MC
paired_lookup::P
dim::D
end
end
12 changes: 5 additions & 7 deletions src/methods/reproject.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,16 +2,15 @@
"""
reproject(obj; crs)
Reproject the lookups (axes) of `obj` to a different crs.
Reproject the lookups of `obj` to a different crs.
This is a lossless operation for the raster data, as only the
lookup values change. This is only possible when the axes of source
and destination projections are aligned: the change is usually from
a [`Regular`](@ref) and an [`Irregular`](@ref) lookup spans.
For converting between projections that are rotated,
skewed or warped in any way, *or* if you want to re-sample the
_data_, use [`resample`](@ref).
skewed or warped in any way, use [`resample`](@ref).
Dimensions without an `AbstractProjected` lookup (such as a `Ti` dimension)
are silently returned without modification.
Expand Down Expand Up @@ -42,9 +41,8 @@ end
"""
reproject(source::GeoFormat, target::GeoFormat, dim::Dimension, val)
`reproject` uses Proj.jl's `Transformation` interface,
but implemented for reprojecting a lookup / axis array,
a single dimension at a time.
`reproject` uses ArchGDAL.reproject, but implemented for a reprojecting
a value array of values, a single dimension at a time.
"""
function reproject(source, target, dim, val)
if source == target
Expand Down Expand Up @@ -81,4 +79,4 @@ end

# Guess the step for arrays
stepof(A::AbstractArray) = (last(A) - first(A)) / (length(A) - 1)
stepof(A::AbstractRange) = step(A)
stepof(A::AbstractRange) = step(A)

0 comments on commit 8ca0e24

Please sign in to comment.