# command to auto-reload modules when they are edited (easier for testing and debugging)
%load_ext autoreload
%autoreload 2

import healpix_geo
import numpy as np
import matplotlib.pyplot as plt
from healpix_analyse.make_rectangle import make_healpix_rectangle_from_lonlat

level = 12
ellipsoid = 'WGS84'

# Coverage box
lon_min, lon_max = 20.0, 40.0
lat_min, lat_max = 40.0, 45.0

cell_ids_2D_array = make_healpix_rectangle_from_lonlat(bbox=(lon_min, lat_min, lon_max, lat_max), 
                                                       level=level, 
                                                       ellipsoid=ellipsoid)

lon, lat = healpix_geo.ring.healpix_to_lonlat(cell_ids_2D_array.flatten(), level, ellipsoid=ellipsoid)

definition = 200
data = np.cos(1024*lon/definition) * np.sin(256*lat/definition) + 1

plt.scatter(lon, lat, c=data, s=1)

<matplotlib.collections.PathCollection at 0x108e7fa30>

cell_ids_2D_array

array([[386418, 386419, 386420, ..., 386514, 386515, 386516],
       [384661, 384662, 384663, ..., 384757, 384758, 384759],
       [382909, 382910, 382911, ..., 383005, 383006, 383007],
       ...,
       [ 62342,  62341,  62340, ...,  62410,  62411,  62412],
       [ 61638,  61637,  61636, ...,  61706,  61707,  61708],
       [ 60938,  60937,  60936, ...,  61005,  61006,  61007]],
      shape=(266, 99), dtype=uint64)

level = 12
ellipsoid = 'WGS84'
indexing_scheme = 'ring' # 'ring' or 'nested'
definition = 20

# Coverage box
lon_min, lon_max = 2.0, 6.0
lat_min, lat_max = 30.0, 32.0

cell_ids, _, _ = healpix_geo.nested.zone_coverage(bbox=(lon_min,lat_min,lon_max,lat_max),
                                                  depth=level, 
                                                  ellipsoid=ellipsoid)

if indexing_scheme == 'ring': 
    cell_ids = to_ring(cell_ids, depth=level)
    lon, lat = healpix_geo.ring.healpix_to_lonlat(cell_ids, level, ellipsoid=ellipsoid)
elif indexing_scheme == 'nested':
    lon, lat = healpix_geo.nested.healpix_to_lonlat(cell_ids, level, ellipsoid=ellipsoid)
else:
    raise NotImplementedError("For now, indexing_scheme must be 'ring' or 'nested'")

data = np.cos(1024*lon/definition) * np.sin(256*lat/definition) + 1

plt.scatter(lon, lat, c=data, s=1)

<matplotlib.collections.PathCollection at 0x303c02770>

nu, cl = powerspectra(cell_ids, level, data, indexing_scheme=indexing_scheme, plot_2D_fft=True)

plt.plot(nu, cl)
plt.yscale('log')
plt.xscale('log')

lats = np.unique(lat)
lons = np.unique(lon)

plt.plot(np.diff(lats), label='lat spacing')
plt.plot(np.diff(lons), label='lon spacing')
plt.legend()
plt.show()