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

try to use healpy.Rotator to apply the rotation #11

Open
wants to merge 5 commits into
base: main
Choose a base branch
from
Open
Changes from 1 commit
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
Prev Previous commit
Next Next commit
try to use healpy.Rotator to apply rot
  • Loading branch information
keewis committed Apr 3, 2024
commit a1956e35fded7f64cee58711d01926f04468e17f
37 changes: 24 additions & 13 deletions xarray_healpy/grid.py
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,9 @@ class HealpyGridInfo:

rot: dict[str, float]

def __post_init__(self):
self.rotator = hp.Rotator(rot=[self.rot[k] for k in ["lon", "lat"]], deg=True)

@property
def nside(self):
return 2**self.level
Expand All @@ -76,21 +79,29 @@ def base_pixels(self, cell_ids):
return np.unique(base_pixel(self.level, cell_ids.data))

def rotate(self, grid, *, direction="rotated"):
if direction == "rotated":
return grid.assign_coords(
{
"latitude": grid["latitude"] - self.rot["lat"],
"longitude": grid["longitude"] - self.rot["lon"],
}
)
elif direction == "global":
return grid.assign_coords(
{
"latitude": grid["latitude"] + self.rot["lat"],
"longitude": grid["longitude"] + self.rot["lon"],
}
def _rotate(lon, lat, inv):
rotated_lon, rotated_lat = self.rotator(
lon.flatten(), lat.flatten(), inv=inv, lonlat=True
)

rotated_lon_ = np.reshape(rotated_lon, lon.shape)
rotated_lat_ = np.reshape(rotated_lat, lat.shape)

return rotated_lon_, rotated_lat_

lon = grid["longitude"]
lat = grid["latitude"]

rotated_lon, rotated_lat = xr.apply_ufunc(
_rotate,
lon,
lat,
kwargs={"inv": direction == "global"},
input_core_dims=[list(lon.dims), list(lat.dims)],
output_core_dims=[list(lon.dims), list(lat.dims)],
)
return grid.assign_coords({"latitude": rotated_lat, "longitude": rotated_lon})

def target_grid(self, source_grid):
if self.rot:
source_grid = self.rotate(source_grid, direction="rotated")
Expand Down