Skip to content
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

Rough methods added get_area_varz and get_rectangle_area_varz #58

Draft
wants to merge 4 commits into
base: develop
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
39 changes: 39 additions & 0 deletions ndsl/grid/gnomonic.py
Original file line number Diff line number Diff line change
Expand Up @@ -385,6 +385,20 @@ def get_area(lon, lat, radius, np):
)


def get_area_varz(lon, lat, radius, height, np):
"""
Given latitude and longitude on cell corners, return the area of each cell.
"""
xyz = lon_lat_to_xyz(lon, lat, np)
lower_left = xyz[(slice(None, -1), slice(None, -1), slice(None, None))]
lower_right = xyz[(slice(1, None), slice(None, -1), slice(None, None))]
upper_left = xyz[(slice(None, -1), slice(1, None), slice(None, None))]
upper_right = xyz[(slice(1, None), slice(1, None), slice(None, None))]
return get_rectangle_area_varz(
lower_left, upper_left, upper_right, lower_right, radius, height, np
)


def set_corner_area_to_triangle_area(
lon, lat, area, tile_partitioner, rank, radius, np
):
Expand Down Expand Up @@ -609,6 +623,31 @@ def get_rectangle_area(p1, p2, p3, p4, radius, np):
return (total_angle - 2 * PI) * radius ** 2


def get_rectangle_area_varz(p1, p2, p3, p4, radius, height, np):
"""
Given four point arrays whose last dimensions are x/y/z in clockwise or
counterclockwise order, return an array of spherical rectangle areas.
NOTE, this is not the exact same order of operations as the Fortran code
This results in some errors in the last digit, but the spherical_angle
is an exact match. The errors in the last digit multipled out by the radius
end up causing relative errors larger than 1e-14, but still wtihin 1e-12.

This method also takes into account the effect of height from the surface of
the Earth on surface area
"""
total_angle = spherical_angle(p2, p3, p1, np)
for (
q1,
q2,
q3,
) in ((p3, p2, p4), (p4, p3, p1), (p1, p4, p2)):
total_angle += spherical_angle(q1, q2, q3, np)

radius = radius + height

return (total_angle - 2 * PI) * radius ** 2


def get_triangle_area(p1, p2, p3, radius, np):
"""
Given three point arrays whose last dimensions are x/y/z, return an array of
Expand Down
2 changes: 1 addition & 1 deletion setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ def local_pkg(name: str, relative_path: str) -> str:

setup(
author="NOAA/NASA",
python_requires=">=3.8",
python_requires=">=3.11.7",
classifiers=[
"Development Status :: 2 - Pre-Alpha",
"Intended Audience :: Developers",
Expand Down
Loading