-
Notifications
You must be signed in to change notification settings - Fork 36
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
optimize cellsize for lon-lat projections #768
base: main
Are you sure you want to change the base?
Conversation
How do we get rid of ArchGDAL? |
I was thinking we could just make a |
But how do we deal with Proj/Wkt, a lot of different text strings can be wgs84 |
The way we already do it - I don't think this line of code needs ArchGDAL
|
Yeah that's ArchGDAL doing that conversion ;) |
I was thrown off by the GeoFormatTypes wrapper, but makes sense. Let's just leave it as it is, then. |
Yes ArchGDAL actually pirates |
Can we take this opportunity to rename |
There is room for significant speed gains for small small grid spacing: #For lat/lon aligned rasters with small grid spacing, this will be much faster and just as accurate:
# load packages
import Rasters: EPSG
using Rasters
using Rasters.Lookups
using DimensionalData
# define meters to lat/lon function
#"""
meters2lonlat_distance(distance_meters, latitude_degrees)
Returns the decimal degree distance along latitude and longitude lines given a distance in
meters and a latitude in decimal degrees.
# Example usage:
#```julia-repl
julia> distance_meters = 1000.0;
julia> latitude_degrees = 45.0;
julia> lat, lon = Altim.meters2lonlat_distance(distance_meters, latitude_degrees)
(0.008997741566866717, 0.012718328120254203)
#```
#"""
function meters2lonlat_distance(distance_meters, latitude_degrees)
# Radius of the Earth in meters
earth_radius = 6371000
# Calculate the angular distance in radians
angular_distance = distance_meters / earth_radius
# Calculate the longitude distance using the Haversine formula
longitude_distance = angular_distance * (180.0 / π) / cosd(latitude_degrees)
latitude_distance = distance_meters / 111139
return latitude_distance, longitude_distance
end
# build an example raster
dX = 0.1
dY = -0.1
lon = X(Projected(166.:dX:168.; sampling=Intervals(Start()), order=ForwardOrdered(), span=Regular(dX), crs=EPSG(4326)))
lat = Y(Projected(-78.0:dY:-80.; sampling=Intervals(Start()), order=ForwardOrdered(), span=Regular(dY ), crs=EPSG(4326)))
ras = Raster(rand(lon, lat))
# given a lat/lon raster with small grid spacing caclculate area
dlon = dims(ras, :X)
dlat = dims(ras, :Y)
lonlat_per_meter = meters2lonlat_distance.(Ref(1), dlat)
dist_lon = step(dlon) ./ getindex.(lonlat_per_meter, 1)
dist_lat = step(dlat) ./ getindex.(lonlat_per_meter, 2)
area = ones(dlon) * DimArray(dist_lon .* dist_lat, dlat)' |
Let's wait for a complete Unitful extension to add the I'm happy to rename but we should |
given the incorrectness of |
Have you compared to the implementation in this PR? For me it is giving similar speeds. I can see the point of your implementation, but using the haversine formula also has downsides, since it doesn't account for the curvature of the earth. I know this isn't a big deal in most cases, but the error increases with the size of gridcells. So after (dis)aggregating a raster, the total area returned by cellarea would be different, which is not ideal. Maybe we should implement it for other projections than lon-lat, but then we might need to do some more geometry. |
Maybe we should already start returning the result in metres though, to avoid another breaking change down the line. If we add a units keyword later then it would default to m^2. |
Yeah result in metres is fine, just not the units kw |
Your implementation knocks the socks off of my proposal (2x faster) and is more accurate so ignore my recommendation. Excited to see this progressing... fast, useful and intuitive... couldn't ask for anything more |
Is this good to go? |
I think so! Can you or @asinghvi17 just nod at the logic in the line of code that reprojects to |
Feel free to ignore the cosd/sind stuff if there's no time, I just discovered that they use an extended-precision deg2rad that might be useful. |
Thanks for reviewing @asinghvi17! I didn't know about cosd and sind so that was really helpful. I think this is really solid now! The only thing I could think of that we could still add is an option to disable the transformation and just return the cartesian area for planar projections. Should I add it? |
I think this would be super useful. If the user understands their projection and just want map ptojected area in map units they could simply pass a kwarg like |
I implemented this now and it works (even without ArchGDAL). I don't know if I like |
Ideally we would say Then if one wants as much precision as possible it's even possible to do |
It seems like now would be the time to do this as this is a breaking change. We could make it easy and just remove kwarg in favor of Linear() and Spherical() types that could be internally defined... then we can expand in future to be more flexible and to eventually integrate GeometryOps... but at a later date? |
It would take about 3 days to get out, but we could technically define a GeometryOpsCore.jl package today and get it registered, Rasters could then depend on that (it should have way less dependencies / load time). I wouldn't want to define them internally because someone using this version of Rasters and a future version of GeometryOps will get quite a bit of incompatibility. |
We have also been talking about singletons like that (in person here, sorry to be exclusive!). But it makes sense to share them with GeometryOps.jl |
ext/RastersArchGDALExt/cellarea.jl
Outdated
_area_from_coords(transform, geom) = _area_from_coords(transform, GI.trait(geom), geom) | ||
function _area_from_coords(transform::AG.CoordTransform, ::GI.LinearRingTrait, ring) | ||
points = map(GI.getpoint(ring)) do p | ||
t = AG.transform!(AG.createpoint(p...), transform) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
createpoint
is super slow, its allocating C objects that need a finalizer and GC for every single point. Can we run the transformation on larger geometries an skip this? Like just GI.convert(LinearRing, ring)
instead?
ext/RastersArchGDALExt/cellarea.jl
Outdated
function _area_from_coords(transform::AG.CoordTransform, ::GI.LinearRingTrait, ring) | ||
points = map(GI.getpoint(ring)) do p | ||
t = AG.transform!(AG.createpoint(p...), transform) | ||
(GI.x(t), GI.y(t)) | ||
end | ||
return _spherical_quadrilateral_area(GI.LinearRing(points)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
function _area_from_coords(transform::AG.CoordTransform, ::GI.LinearRingTrait, ring) | |
points = map(GI.getpoint(ring)) do p | |
t = AG.transform!(AG.createpoint(p...), transform) | |
(GI.x(t), GI.y(t)) | |
end | |
return _spherical_quadrilateral_area(GI.LinearRing(points)) | |
function _area_from_coords(transform::AG.CoordTransform, trait::GI.AbstractCurveTrait, ring) | |
t = AG.transform!(GI.convert(AG.geointerface_geomtype(trait), ring), transform) | |
return _spherical_quadrilateral_area(GI.convert(GI.geointerface_geomtype(trait), t)) |
Something like this should be better as GI skips around ArchGDAL point creation and just works on the internal vectors of points in the LinearRing
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Guessing this will be 3x faster
It also means a lot less calls to transform!
as well as less allocations
(We could also just add a Proj dep and skip around GDAL completely, Proj is super fast just on points...)
ext/RastersArchGDALExt/cellarea.jl
Outdated
GI.LinearRing([ | ||
(xb[1], yb[1]), | ||
(xb[2], yb[1]), | ||
(xb[2], yb[2]), | ||
(xb[1], yb[2]), | ||
(xb[1], yb[1]) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This might give some more micro optimisation, 4 things is usually better than 5 things with computers...
GI.LinearRing([ | |
(xb[1], yb[1]), | |
(xb[2], yb[1]), | |
(xb[2], yb[2]), | |
(xb[1], yb[2]), | |
(xb[1], yb[1]) | |
GI.LineSting([ | |
(xb[1], yb[1]), | |
(xb[2], yb[1]), | |
(xb[2], yb[2]), | |
(xb[1], yb[2]), |
ext/RastersArchGDALExt/cellarea.jl
Outdated
|
||
function _spherical_quadrilateral_area(ring) | ||
ps = GI.getpoint(ring) | ||
(p1, p2, p3, p4) = _lonlat_to_sphericalpoint.((ps[1], ps[2], ps[3], ps[4])) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The last point isn't actually used here so my LineString
approach below does make sense
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #768 +/- ##
==========================================
+ Coverage 82.32% 82.61% +0.28%
==========================================
Files 60 62 +2
Lines 4357 4566 +209
==========================================
+ Hits 3587 3772 +185
- Misses 770 794 +24 ☔ View full report in Codecov by Sentry. |
Happy to report that this solves #764 in web mercator as well! |
Solves #764 and makes cellsize calculation much faster for lon-lat projections.
I'm considering getting rid of the ArchGDAL requirement for these projections as well.
Any ideas for further improvements @rafaqz @alex-s-gardner?