From e741d2af27b7aa29edfa1fd7f00b513150578c8d Mon Sep 17 00:00:00 2001 From: tiemvanderdeure Date: Wed, 30 Aug 2023 18:08:14 +0200 Subject: [PATCH 1/4] add cellsize function --- ext/RastersArchGDALExt/RastersArchGDALExt.jl | 3 +- ext/RastersArchGDALExt/cellsize.jl | 109 +++++++++++++++++++ src/Rasters.jl | 2 +- src/extensions.jl | 29 ++--- test/cellsize.jl | 25 +++++ 5 files changed, 148 insertions(+), 20 deletions(-) create mode 100644 ext/RastersArchGDALExt/cellsize.jl create mode 100644 test/cellsize.jl diff --git a/ext/RastersArchGDALExt/RastersArchGDALExt.jl b/ext/RastersArchGDALExt/RastersArchGDALExt.jl index 32e568e5e..63dce90e2 100644 --- a/ext/RastersArchGDALExt/RastersArchGDALExt.jl +++ b/ext/RastersArchGDALExt/RastersArchGDALExt.jl @@ -20,7 +20,7 @@ using Rasters: GDALsource, AbstractProjected, RasterStackOrArray, FileArray, RES_KEYWORD, SIZE_KEYWORD, CRS_KEYWORD, FILENAME_KEYWORD, SUFFIX_KEYWORD, EXPERIMENTAL, GDAL_EMPTY_TRANSFORM, GDAL_TOPLEFT_X, GDAL_WE_RES, GDAL_ROT1, GDAL_TOPLEFT_Y, GDAL_ROT2, GDAL_NS_RES -import Rasters: reproject, resample, warp +import Rasters: reproject, resample, warp, cellsize const RA = Rasters const DD = DimensionalData @@ -28,6 +28,7 @@ const DA = DiskArrays const GI = GeoInterface const LA = LookupArrays +include("cellsize.jl") include("gdal_source.jl") include("reproject.jl") include("resample.jl") diff --git a/ext/RastersArchGDALExt/cellsize.jl b/ext/RastersArchGDALExt/cellsize.jl new file mode 100644 index 000000000..51b855012 --- /dev/null +++ b/ext/RastersArchGDALExt/cellsize.jl @@ -0,0 +1,109 @@ +""" + cellsize(x) + +Gives the approximate size of each cell in square km. +This function works for any projection, but uses an algorithm for spherical polygons, and thus gives the approximates the true size to about 0.1%, depending on latitude. + +# Arguments + +- `x`: A `Raster` or a `Tuple` of `X` and `Y` dimensions. + +## Example + +```jldoctest +using Rasters, Rasters.LookupArrays +dimz = X(Projected(90.0:0.1:120; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326))), + Y(Projected(0.0:0.1:50; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326))) + +cs = cellsize(dimz) + +``` +$EXPERIMENTAL +""" + +function _great_circle_bearing(lon1::AbstractFloat, lat1::AbstractFloat, lon2::AbstractFloat, lat2::AbstractFloat) + dLong = lon1 - lon2 + + s = cos(lat2)*sin(dLong) + c = cos(lat1)*sin(lat2) - sin(lat1)*cos(lat2)*cos(dLong) + + return atan(s, c) +end + +## Get the area of a LinearRing with coordinates in radians +# Using Gidard's theorem +function _area_from_rads(ring; R = 6371.0088) + n = GI.npoint(ring) + area = -(n-3)*pi + + prevpoint = GI.getpoint(ring, n-1) + point = GI.getpoint(ring, 1) + + for i in 2:n + nextpoint = GI.getpoint(ring, i) + + beta1 = _great_circle_bearing(GI.x(point), GI.y(point), GI.x(prevpoint), GI.y(prevpoint)) + beta2 = _great_circle_bearing(GI.x(point), GI.y(point), GI.x(nextpoint), GI.y(nextpoint)) + angle = acos(cos(-beta1)*cos(-beta2) + sin(-beta1)*sin(-beta2)) + area += angle + + prevpoint = point + point = nextpoint + end + + return area*R^2 +end + +_area_from_coords(transform, geom) = _area_from_coords(transform, GI.trait(geom), geom) +_area_from_coords(geom) = _area_from_coords(GI.trait(geom), geom) + +function _area_from_coords(::GI.LinearRingTrait, ring) # no ArchGDAL, assumes degrees + points = map(GI.getpoint(ring)) do p + (deg2rad(GI.x(p)), deg2rad(GI.y(p))) + end + + return _area_from_rads(GI.LinearRing(points)) +end + +function _area_from_coords(transform::ArchGDAL.CoordTransform, ::GI.LinearRingTrait, ring) + points = map(GI.getpoint(ring)) do p + t = ArchGDAL.transform!(ArchGDAL.createpoint(p...), transform) + (deg2rad(GI.x(t)), deg2rad(GI.y(t))) + end + return _area_from_rads(GI.LinearRing(points)) +end + +function cellsize(x::Tuple{X, Y}) + xbnds, ybnds = DimensionalData.intervalbounds(x) + + if convert(CoordSys, crs(x)) == CoordSys("Earth Projection 1, 104") # check if need to reproject + areas = [_area_from_coords( + GI.LinearRing([ + (xb[1], yb[1]), + (xb[2], yb[1]), + (xb[2], yb[2]), + (xb[1], yb[2]), + (xb[1], yb[1]) + ])) + for xb in xbnds, yb in ybnds] + else + areas = ArchGDAL.crs2transform(crs(x), EPSG(4326)) do transform + [_area_from_coords( + transform, + GI.LinearRing([ + (xb[1], yb[1]), + (xb[2], yb[1]), + (xb[2], yb[2]), + (xb[1], yb[2]), + (xb[1], yb[1]) + ])) + for xb in xbnds, yb in ybnds] + end + end + + return Raster(areas, dims = x) +end + +function cellsize(x::Raster) + cellsize(dims(x, (X, Y))) +end diff --git a/src/Rasters.jl b/src/Rasters.jl index 0af7c82ab..a982badb0 100644 --- a/src/Rasters.jl +++ b/src/Rasters.jl @@ -60,7 +60,7 @@ export missingval, boolmask, missingmask, replace_missing, replace_missing!, aggregate, aggregate!, disaggregate, disaggregate!, mask, mask!, resample, warp, zonal, crop, extend, trim, slice, combine, points, classify, classify!, mosaic, mosaic!, extract, rasterize, rasterize!, - coverage, coverage!, setcrs, setmappedcrs, smapseries + coverage, coverage!, setcrs, setmappedcrs, smapseries, cellsize export crs, mappedcrs, mappedindex, mappedbounds, projectedindex, projectedbounds export reproject, convertlookup diff --git a/src/extensions.jl b/src/extensions.jl index 27efd4c91..6bdd4c008 100644 --- a/src/extensions.jl +++ b/src/extensions.jl @@ -1,29 +1,22 @@ # extensions - -# stubs that need ArchGDAL -function resample(args...; kw...) +function throw_extention_error(package::String, extention::Symbol) @static if isdefined(Base, :get_extension) # julia > 1.9 - if isnothing(Base.get_extension(Rasters, :RastersArchGDALExt)) - throw(BackendException("ArchGDAL")) - else - throw(MethodError(resample, args)) - end + if isnothing(Base.get_extension(Rasters, extention)) + throw(BackendException(package)) else - throw(BackendException("ArchGDAL")) + throw(MethodError(resample, args)) end -end -function warp(args...; kw...) - @static if isdefined(Base, :get_extension) # julia > 1.9 - if isnothing(Base.get_extension(Rasters, :RastersArchGDALExt)) - throw(BackendException("ArchGDAL")) - else - throw(MethodError(warp, args)) - end else - throw(BackendException("ArchGDAL")) + throw(BackendException(package)) end end + +# stubs that need ArchGDAL +resample(args...; kw...) = throw_extention_error("ArchGDAL", :RastersArchGDALExt) +warp(args...; kw...) = throw_extention_error("ArchGDAL", :RastersArchGDALExt) +cellsize(args...; kw...) = throw_extention_error("ArchGDAL", :RastersArchGDALExt) + # Other shared stubs function layerkeys end function smapseries end diff --git a/test/cellsize.jl b/test/cellsize.jl new file mode 100644 index 000000000..f5bea825c --- /dev/null +++ b/test/cellsize.jl @@ -0,0 +1,25 @@ +using Rasters, Rasters.LookupArrays, ArchGDAL +using Test + +include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl")) + +@testset "cellsize" begin + dimz = X(Projected(90.0:0.1:99.9; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326))), + Y(Projected(0.0:0.1:89.9; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326))) + + dimz_25832 = X(Projected(0.0:100:10000.0; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(100), crs=EPSG(25832))), + Y(Projected(0.0:100:10000.0; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(100), crs=EPSG(25832))) + + cs = cellsize(dimz) + cs2 = cellsize(dimz_25832) + + # Check the output is a raster + @test cs isa Raster + @test cs2 isa Raster + # Test that the total area matches the expected area (1/72th of the Earth surface) + @test sum(cs) ≈ 510.1e6/72 rtol = 0.01 + # Test all areas are about 0.01 km2 + @test maximum(cs2) ≈ 0.01 rtol = 0.01 + @test minimum(cs2) ≈ 0.01 rtol = 0.01 + +end From ce71ee5339ce4111068da3dd485a9a7c27ed40b2 Mon Sep 17 00:00:00 2001 From: tiemvanderdeure Date: Thu, 31 Aug 2023 08:20:06 +0200 Subject: [PATCH 2/4] fix typo --- src/extensions.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/extensions.jl b/src/extensions.jl index 6bdd4c008..fae40342f 100644 --- a/src/extensions.jl +++ b/src/extensions.jl @@ -1,7 +1,7 @@ # extensions -function throw_extention_error(package::String, extention::Symbol) +function throw_extension_error(package::String, extension::Symbol) @static if isdefined(Base, :get_extension) # julia > 1.9 - if isnothing(Base.get_extension(Rasters, extention)) + if isnothing(Base.get_extension(Rasters, extension)) throw(BackendException(package)) else throw(MethodError(resample, args)) @@ -13,9 +13,9 @@ end # stubs that need ArchGDAL -resample(args...; kw...) = throw_extention_error("ArchGDAL", :RastersArchGDALExt) -warp(args...; kw...) = throw_extention_error("ArchGDAL", :RastersArchGDALExt) -cellsize(args...; kw...) = throw_extention_error("ArchGDAL", :RastersArchGDALExt) +resample(args...; kw...) = throw_extension_error("ArchGDAL", :RastersArchGDALExt) +warp(args...; kw...) = throw_extension_error("ArchGDAL", :RastersArchGDALExt) +cellsize(args...; kw...) = throw_extension_error("ArchGDAL", :RastersArchGDALExt) # Other shared stubs function layerkeys end From 4729abf0bdb2e7430ea339aca1775da41cb6ecad Mon Sep 17 00:00:00 2001 From: tiemvanderdeure Date: Thu, 31 Aug 2023 10:18:29 +0200 Subject: [PATCH 3/4] fix extension error function --- src/extensions.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/extensions.jl b/src/extensions.jl index fae40342f..0fdf95afe 100644 --- a/src/extensions.jl +++ b/src/extensions.jl @@ -1,10 +1,10 @@ # extensions -function throw_extension_error(package::String, extension::Symbol) +function throw_extension_error(f::Function, package::String, extension::Symbol, args) @static if isdefined(Base, :get_extension) # julia > 1.9 if isnothing(Base.get_extension(Rasters, extension)) throw(BackendException(package)) else - throw(MethodError(resample, args)) + throw(MethodError(f, args)) end else throw(BackendException(package)) @@ -13,9 +13,9 @@ end # stubs that need ArchGDAL -resample(args...; kw...) = throw_extension_error("ArchGDAL", :RastersArchGDALExt) -warp(args...; kw...) = throw_extension_error("ArchGDAL", :RastersArchGDALExt) -cellsize(args...; kw...) = throw_extension_error("ArchGDAL", :RastersArchGDALExt) +resample(args...; kw...) = throw_extension_error(resample, "ArchGDAL", :RastersArchGDALExt, args) +warp(args...; kw...) = throw_extension_error(warp, "ArchGDAL", :RastersArchGDALExt, args) +cellsize(args...; kw...) = throw_extension_error(cellsize, "ArchGDAL", :RastersArchGDALExt, args) # Other shared stubs function layerkeys end From 0f41810f7e09f04e7e1d18831d1b4905dd944769 Mon Sep 17 00:00:00 2001 From: tiemvanderdeure Date: Thu, 31 Aug 2023 10:43:06 +0200 Subject: [PATCH 4/4] minor changes to docstring and formatting --- ext/RastersArchGDALExt/cellsize.jl | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/ext/RastersArchGDALExt/cellsize.jl b/ext/RastersArchGDALExt/cellsize.jl index 51b855012..833008d3c 100644 --- a/ext/RastersArchGDALExt/cellsize.jl +++ b/ext/RastersArchGDALExt/cellsize.jl @@ -2,7 +2,7 @@ cellsize(x) Gives the approximate size of each cell in square km. -This function works for any projection, but uses an algorithm for spherical polygons, and thus gives the approximates the true size to about 0.1%, depending on latitude. +This function works for any projection, using an algorithm for polygons on a sphere. It approximates the true size to about 0.1%, depending on latitude. # Arguments @@ -56,7 +56,6 @@ end _area_from_coords(transform, geom) = _area_from_coords(transform, GI.trait(geom), geom) _area_from_coords(geom) = _area_from_coords(GI.trait(geom), geom) - function _area_from_coords(::GI.LinearRingTrait, ring) # no ArchGDAL, assumes degrees points = map(GI.getpoint(ring)) do p (deg2rad(GI.x(p)), deg2rad(GI.y(p))) @@ -64,7 +63,6 @@ function _area_from_coords(::GI.LinearRingTrait, ring) # no ArchGDAL, assumes de return _area_from_rads(GI.LinearRing(points)) end - function _area_from_coords(transform::ArchGDAL.CoordTransform, ::GI.LinearRingTrait, ring) points = map(GI.getpoint(ring)) do p t = ArchGDAL.transform!(ArchGDAL.createpoint(p...), transform) @@ -73,10 +71,9 @@ function _area_from_coords(transform::ArchGDAL.CoordTransform, ::GI.LinearRingTr return _area_from_rads(GI.LinearRing(points)) end -function cellsize(x::Tuple{X, Y}) - xbnds, ybnds = DimensionalData.intervalbounds(x) - - if convert(CoordSys, crs(x)) == CoordSys("Earth Projection 1, 104") # check if need to reproject +function cellsize(dims::Tuple{X, Y}) + xbnds, ybnds = DimensionalData.intervalbounds(dims) + if convert(CoordSys, crs(dims)) == CoordSys("Earth Projection 1, 104") # check if need to reproject areas = [_area_from_coords( GI.LinearRing([ (xb[1], yb[1]), @@ -87,7 +84,7 @@ function cellsize(x::Tuple{X, Y}) ])) for xb in xbnds, yb in ybnds] else - areas = ArchGDAL.crs2transform(crs(x), EPSG(4326)) do transform + areas = ArchGDAL.crs2transform(crs(dims), EPSG(4326)) do transform [_area_from_coords( transform, GI.LinearRing([ @@ -101,7 +98,7 @@ function cellsize(x::Tuple{X, Y}) end end - return Raster(areas, dims = x) + return Raster(areas, dims) end function cellsize(x::Raster)