''' A gridded sphere.
.. codeauthor:: Yoshifumi Futaana
Sphere gridded into the longitude-latitude system is represented by :class:`SimpleGridSphere`.
Similar to :mod:`tessellatedsphere`, which does not support any container functionality.
Similar to :mod:`polarframe` which support 3D grid and volume.
'''
import logging
_logger = logging.getLogger(__name__)
import pickle as pickle
import bisect
import numpy as np
[docs]class SimpleGridSphere:
''' A simple gridding using longitude latitude system.
A simple container in (pseudo) 2-D gridded map.
Data on geometric map can be saved in it.
A grid sphere object is generated by calling
:class:`SimpleGridSphere` class as following.
>>> sg = SimpleGridSphere(nlon=360, nlat=180)
*Longitude/Latitude pair and index*
You can specify the number of grids both longitude and latitude directions.
Any data on a map can be assigned to the sphere surface.
To specify the coordinates, a pair of longitude and latitude can be used.
Internally, the pair is converted to a single, unique index for each grid.
The index is also usedto specify the coordinate.
The conversion from longitude and latitude to the index can be done via :meth:`lonlat2index`.
>>> sg.lonlat2index(35.78, 28.99)
42515
This means the surface grid containing the coordinates of (lon=35.78, lat=28.99) has index of 42515.
Oppositely, the index to the center of the longitude, latitude is also available (:meth:`index2lonlat`).
For example, the lon and lat of surface element (center location) with id=2500 is 340.5 and -83.5, as follows.
>>> sg.index2lonlat(2500)
(340.5, -83.5)
Similar conversion from index to xyz vector is also available (:meth:`index2xyz`).
*Value on a map*
Any data (usually scalar) can be stored by :meth:`append_value_lonlat` or :meth:`append_value` methods.
The formar takes the pair of longitude and latitude to specify the coordinate,
while the latter takes the index.
>>> sg.append_value_lonlat(25.77, -50.27, -275.25) # Put the data -275.25 to lon=25.77 and lat=-50.27
>>> sg.append_value(2500, 275.25)
*Access to the stored data*
The direct access to the data by index is possible.
>>> print(sg[2500]) # Access by index
[275.25]
>>> print(sg[sg.lonlat2index(25.77, -50.27)]) # Access using lon/lat pair
[-275.25]
Indeed, any type of values can be stored; while it is recommended to use
"numerical" data to store, to work with various methods.
>>> sg.append_value(3500, 'abcdefg')
>>> sg.append_value(3500, 250.5)
>>> print(sg[3500])
['abcdefg', 250.5]
The :meth:`get_average` and :meth:`get_counts` methods will return the
2-D array of the data. For usability, the data is already re-shaped into 2-D array.
Note that the order of the index is ``[latitude][longitude]``.
>>> counts = sg.get_counts()
>>> counts.shape
(180, 360)
>>> print(counts.sum())
4
Indeed you have set a string into the data, thus, the element at index 3500
is set to np.nan when you get the average with warning message.
>>> avg = sg.get_average()
>>> print(avg.shape)
(180, 360)
Manupulations can be available :meth:`get_max`, :meth:`get_min` and :meth:`get_median`
in addition to :meth:`get_average`.
>>> lon, lat = sg.index2lonlat(0)
>>> sg.append_value_lonlat(lon, lat, 30.)
>>> sg.append_value_lonlat(lon, lat, 31.)
>>> sg.append_value_lonlat(lon, lat, 32.)
>>> sg.append_value_lonlat(lon, lat, 33.)
>>> sg.append_value_lonlat(lon, lat, 34.)
>>> sg.append_value_lonlat(lon, lat, 35.)
>>> sg.append_value_lonlat(lon, lat, 43.)
Average should be 34, median to be 33, min to be 30 and max to be 43.
>>> avg = sg.get_average()
>>> print(avg[0, 0])
34.0
>>> medi = sg.get_median()
>>> print(medi[0, 0])
33.0
>>> maxx = sg.get_max()
>>> print(maxx[0, 0])
43.0
>>> minn = sg.get_min()
>>> print(minn[0, 0])
30.0
*Plotting*
For showing by pcolor in the matplotlib, you can just get the grid values by
:meth:`get_cgrid` or :meth:`get_bgrid`.
The first one returns the central location, but the second one returns the boundary locations.
>>> lonc, latc = sg.get_cgrid()
>>> print(lonc.shape)
(180, 360)
>>> lonb, latb = sg.get_bgrid()
>>> print(lonb.shape)
(181, 361)
To plot using pseudo-color is useful using the results of :meth:`get_bgrid`.
.. code-block:: py
plt.pcolor(lonb, latb, avg)
*Pickling*
Pickling can be done.
>>> import pickle as pickle
>>> pickled = pickle.dumps(sg)
>>> sg2 = pickle.loads(pickled)
>>> sg is sg2
False
>>> sg.datalist == sg2.datalist
True
>>> sg.datalist is sg2.datalist
False
'''
_logger = logging.getLogger(__name__ + '.SimpleGridSphere')
def __init__(self, nlon=360, nlat=180):
lonlist = np.linspace(0., 360., nlon+1)
latlist = np.linspace(-90., 90., nlat+1)
lon0list = lonlist[:-1]
lon1list = lonlist[1:]
lat0list = latlist[:-1]
lat1list = latlist[1:]
lonclist = (lon0list + lon1list) / 2.
latclist = (lat0list + lat1list) / 2.
lon0, lat0 = np.meshgrid(lon0list, lat0list) # Shape lons is (nlat+1, nlon+1)
lon1, lat1 = np.meshgrid(lon1list, lat1list) # Shape lons is (nlat+1, nlon+1)
lonc, latc = np.meshgrid(lonclist, latclist) # Shape lons is (nlat+1, nlon+1)
self.lonc = lonc.flatten()
''' Longitude list (flattened into 1D).
'''
self.latc = latc.flatten()
''' Latitude list (flattened into 1D).
'''
self.vectorc = np.array([np.cos(self.latc * np.pi / 180.) * np.cos(self.lonc * np.pi / 180.),
np.cos(self.latc * np.pi / 180.) * np.sin(self.lonc * np.pi / 180.),
np.sin(self.latc * np.pi / 180.)])
''' Central position in vector form. ``np.array`` with (3, N).
'''
self.lon0 = lon0.flatten()
''' Longitude bound list'''
self.lon1 = lon1.flatten()
''' Longitude bound list'''
self.lat0 = lat0.flatten()
''' Latitude bound list'''
self.lat1 = lat1.flatten()
''' Latitude bound list'''
self.nlat = nlat
''' Number of latitude '''
self.nlon = nlon
''' Number of longitude'''
self.datalist = [[] for i in self.lonc]
''' Datalist. Array with size of ``N``. Each element is also array'''
self.angular_distance_matrix = None
''' A test feature of angular distance matrix.
Not recommended to use. See :meth:`make_angular_distance_matrix`'''
def __len__(self):
return len(self.datalist)
[docs] def index2lonlat(self, index):
''' Return the (lonc, latc) for the given index.
>>> gsph = SimpleGridSphere(nlon=18, nlat=9)
>>> print(gsph.index2lonlat(120))
(250.0, 40.0)
>>> print(gsph.index2lonlat([120, 121, 122])) # doctest: +NORMALIZE_WHITESPACE
(array([250., 270., 290.]), array([40., 40., 40.]))
'''
return (self.lonc[index], self.latc[index])
[docs] def index2xyz(self, index):
return self.vectorc[:, index]
[docs] def lonlat2index(self, lon, lat):
''' Return the index corresponding to (lon, lat), in degrees.
>>> m = SimpleGridSphere()
>>> print(m.lonlat2index(0, -90))
0
>>> print(m.lonlat2index(1, -90))
1
>>> print(m.lonlat2index(1.5, -90))
1
>>> print(m.lonlat2index(360, -90))
0
>>> print(m.lonlat2index(360, -89))
360
>>> print(m.lonlat2index(-360, 89))
64440
>>> print(m.lonlat2index(0, 88.9))
64080
>>> print(m.lonlat2index(0, 90))
64440
'''
if not (-90 <= lat <= 90):
raise ValueError('lat (%f) out of range (-90<=lat<=90)' % lat)
# if lat == 90:
# lat = self.latc.max() # 90 will give error in the following. Instead, replace it to the maximum center latitude.
if not np.isfinite(lon):
raise ValueError('lon {} should be finate value.'.format(lon))
while lon < 0:
lon += 360
while lon >= 360:
lon -= 360
# lonlist = np.linspace(0., 360., self.nlon+1)
# latlist = np.linspace(-90., 90., self.nlat+1)
ilon = bisect.bisect(self.lon0[:self.nlon], lon) - 1
ilat = bisect.bisect(self.lat0[::self.nlon], lat) - 1
if ilon < 0 or ilon >= self.nlon:
raise ValueError('Something wrong in longitude index. %f' % lon)
if ilat < 0 or ilat >= self.nlat:
raise ValueError('Something wrong in latgitude index. %f' % lat)
idx = ilon + ilat * self.nlon
return idx
[docs] def lonlats2indexes(self, lon_list, lat_list):
lonlist = np.array(lon_list)
latlist = np.array(lat_list)
if latlist.max() > 90 or latlist.min() < -90:
raise ValueError('latlist out of range [-90, 90]. %f %f.'
% (latlist.max(), latlist.min()))
while lonlist.min() < 0:
lonlist = np.where(lonlist < 0, lonlist + 360, lonlist)
while lonlist.max() >= 360:
lonlist = np.where(lonlist >= 360, lonlist - 360, lonlist)
lon0 = self.lon0[:self.nlon]
lat0 = self.lat0[::self.nlon]
ilonlist = np.array([bisect.bisect(lon0, lon) -1 for lon in lonlist])
ilatlist = np.array([bisect.bisect(lat0, lat) -1 for lat in latlist])
if ilonlist.min() < 0 or ilonlist.max() >= self.nlon:
raise ValueError('Something wrong in longitude index.')
if ilatlist.min() < 0 or ilatlist.max() >= self.nlat:
raise ValueError('Something wrong in latgitude index.')
idxlist = ilonlist + ilatlist * self.nlon
return idxlist
[docs] def get_cgrid(self):
''' Return the center grid as a pair of 2-D np.array.
:returns: (``lonarray``, ``latarray``) where ``lonarray`` (``latarray``) is the center
longitudal (latitudal) grid in (``nlat``, ``nlon``) shape.
'''
return (self.lonc.reshape([self.nlat, self.nlon]), self.latc.reshape([self.nlat, self.nlon]))
[docs] def get_bgrid(self, delta=0):
''' Return the boundary grid as a pair of 2-D np.array.
:keyword delta: The edges of longitude and latitude will be shrinkied with this value.
If the default value 0.001 is used, the longitude range is from 0.001 to 359.999
and the latitude range is from -89.999 to 89.999.
This is a bad trick, but s very useful for plotting.
:returns: (``lonarray``, ``latarray``) where ``lonarray`` (``latarray``) is the boundary
longitudal (latitudal) grid in (``nlat+1``, ``nlon+1``) shape.
'''
# Simpler to re-generate them...
lonlist = np.linspace(0. + delta, 360. - delta, self.nlon + 1)
latlist = np.linspace(-90. + delta, 90. - delta, self.nlat + 1)
return np.meshgrid(lonlist, latlist)
[docs] def get_average(self):
''' Get average of each grid.
:return: ``np.ma.masked_array`` with shape of (nlat, nlon)
:rtype: ``np.ma.masked_array``
It works correctly only when all the data is ``float`` (or equivalent).
Otherwise, masked.
>>> gs = SimpleGridSphere()
>>> gs.append_value_lonlat(150, 65, 25.3)
>>> gs.append_value_lonlat(358, 27.5, 883)
>>> gs.append_value_lonlat(25, 90, -35)
>>> lons, lats = gs.get_bgrid()
>>> ave = gs.get_average()
You can use pcolor with::
plt.pcolor(lons, lats, ave)
'''
return self.__get_something(np.mean)
[docs] def get_average2(self):
''' Return the averaged value of the data in each grid.
It is equivalent to :meth:`get_average`, but more primitive implementaiton.
Recommended to use :meth:`get_average` for usual purpose.
'''
avelist = np.empty([len(self.datalist)])
avelist.fill(np.nan)
for index, data in enumerate(self.datalist):
dataarray = np.array(data)
try:
if dataarray.size == 0:
continue
average = dataarray.mean()
except TypeError as e:
self._logger.warning('Non-float value had been added at index %d. (%s)' % (index, str(e)))
continue
avelist[index] = average
avelist.shape = [self.nlat, self.nlon]
return np.ma.masked_invalid(avelist)
def __get_something(self, evalfunc):
''' A function that returns something
'''
somethinglist = np.empty([len(self.datalist)])
somethinglist.fill(np.nan)
for index, data in enumerate(self.datalist):
dataarray = np.array(data)
try:
if dataarray.size == 0:
continue
something = evalfunc(dataarray)
except TypeError as e:
self._logger.warning('No calculation can be done.')
continue
somethinglist[index] = something
somethinglist.shape = [self.nlat, self.nlon]
return np.ma.masked_invalid(somethinglist)
[docs] def get_max(self):
return self.__get_something(np.max)
[docs] def get_min(self):
return self.__get_something(np.min)
[docs] def get_counts(self):
''' Return the number of data stored in the grid
Similar to :meth:`get_average`, while usual np.array is returned.
'''
retlist = np.zeros([len(self.datalist)], dtype=int)
for index, data in enumerate(self.datalist):
retlist[index] = len(data)
retlist.shape = [self.nlat, self.nlon]
return retlist
[docs] def append_value(self, index, value):
''' Append the given value to the given index.
'''
self.datalist[index].append(value)
[docs] def append_values(self, indexlist, valuelist):
for i, v in zip(indexlist, valuelist):
self.append_value(i, v)
[docs] def append_value_lonlat(self, lon, lat, value):
''' Append the given value to the given (lon, lat).
:param lon: Longitude (degrees)
:param lat: Latitude (degrees)
:param value: Value
'''
idx = self.lonlat2index(lon, lat)
self.append_value(idx, value)
[docs] def append_values_lonlat(self, lonlist, latlist, valuelist):
idxlist = self.lonlats2indexes(lonlist, latlist)
self.append_values(idxlist, valuelist)
def __getitem__(self, index):
''' Index-based access getter support.
'''
return self.datalist[index]
def __add__(self, other):
''' Adding two SimpleGridSphere.
>>> gs1 = SimpleGridSphere(nlon=4, nlat=6)
>>> gs2 = SimpleGridSphere(nlon=4, nlat=6)
>>> gs1.append_value(5, 15.)
>>> gs1.append_value(5, 16)
>>> gs1.append_value(5, 17)
>>> print(gs1.get_counts().sum())
3
>>> gs2.append_value(5, 33)
>>> gs2.append_value(5, 34)
>>> print(gs2.get_counts().sum())
2
>>> gs = gs1 + gs2
>>> print(len(gs))
24
>>> print(gs.get_counts().sum())
5
>>> print(gs[5])
[15.0, 16, 17, 33, 34]
>>> print(gs.get_average()[1, 1]) #(1, 1) is flattened index of 5
23.0
'''
if self.nlon != other.nlon:
raise ValueError('Longitudal separation is not the same %d&%d'
% (self.nlon, other.nlon))
elif self.nlat != other.nlat:
raise ValueError('Latitude separation is not the same %d&%d'
% (self.nlat, other.nlat))
added = SimpleGridSphere(nlon=self.nlon, nlat=self.nlat)
# for i in range(len(added)):
# for s in self.datalist[i]:
# added.datalist[i].append(s)
# for s in other.datalist[i]:
# added.datalist[i].append(s)
for i in range(len(added)):
added.datalist[i] = self.datalist[i] + other.datalist[i]
return added
[docs] def clear(self):
''' Clear the data. All the list in the data will be replaced by ``[]``.
'''
self.datalist = [[] for i in self.lonc]
self.make_angular_distance_matrix(clear_cache=True)
[docs] def make_angular_distance_matrix(self, sqrt=True, clear_cache=False):
''' Create a distance matrix.
The distance matrix is the matrix that saves the distance between the elements.
It can be used, for example, to calculate the spacecraft visibility from the surface of body.
.. warning::
The method will consume time and memory quite a lot.
This functionality is memory hungry. Indeed, it consumes (``nlon`` * ``nlat``) ** 2.
This means that, as default gridding (360, 180) will need 4.2e9 data, which may consume
more than 10 GB memory!!
However, for smaller gridding, e.g. (36, 18) every 10 degrees, may reduces 4.2e5 data,
which may consume only a few MB.
This functionality also time consuming.
On iMac in my office it takes ~1 sec for nlon*nlat = 144.
The algorithm is N(O^2), thus, (nlon=36, nlat=18) will take ~25 sec.
Generated matrix will have (N, N) elements, where N is nlon * nlat.
``matrix[i, j]`` have the minimum distance between the grid indexed i and j for i < j.
``matrix[i, j]`` have the maximum distance between the grid indexed i and j for i > j.
``matrix[i, j]`` have zero for i == j.
>>> gs = SimpleGridSphere(nlon=6, nlat=3)
>>> gs.make_angular_distance_matrix()
'''
# from pudb import set_trace; set_trace()
if clear_cache:
self.angular_distance_matrix = None
return
n = len(self)
adm = np.ones([n, n])
a2v = lambda lonr, latr: np.array([np.cos(lonr) * np.cos(latr), np.sin(lonr) * np.cos(latr), np.sin(latr)])
# First, calculate the vector for each point.
vectors = [] # (n, 4, 3) array
for i in range(n):
lon0 = self.lon0[i] * np.pi / 180.
lon1 = self.lon1[i] * np.pi / 180.
lat0 = self.lat0[i] * np.pi / 180.
lat1 = self.lat1[i] * np.pi / 180.
v0 = a2v(lon0, lat0)
v1 = a2v(lon0, lat1)
v2 = a2v(lon1, lat0)
v3 = a2v(lon1, lat1)
vectors.append([v0, v1, v2, v3])
vectors = np.array(vectors)
for i in range(n):
# adm[i, i] = 1 # Not needed. Initilized by 1.
vec_i = vectors[i, :, :] # (4, 3) shape
for j in range(i+1, n): # i < j, always.
vec_j = vectors[j, :, :] # (4, 3) shape
angle16 = np.empty([16]) # Set 16 pair of angles
#
angle16[0:4] = vec_i[0, :].dot(vec_j.T)
angle16[4:8] = vec_i[1, :].dot(vec_j.T)
angle16[8:12] = vec_i[2, :].dot(vec_j.T)
angle16[12:16] = vec_i[3, :].dot(vec_j.T)
# Equivalent to followings.
# # To be more simplifid(?)
# angle16[0] = vec_i[0, :].dot(vec_j[0, :])
# angle16[1] = vec_i[0, :].dot(vec_j[1, :])
# angle16[2] = vec_i[0, :].dot(vec_j[2, :])
# angle16[3] = vec_i[0, :].dot(vec_j[3, :])
# angle16[4] = vec_i[1, :].dot(vec_j[0, :])
# angle16[5] = vec_i[1, :].dot(vec_j[1, :])
# angle16[6] = vec_i[1, :].dot(vec_j[2, :])
# angle16[7] = vec_i[1, :].dot(vec_j[3, :])
# angle16[8] = vec_i[2, :].dot(vec_j[0, :])
# angle16[9] = vec_i[2, :].dot(vec_j[1, :])
# angle16[10] = vec_i[2, :].dot(vec_j[2, :])
# angle16[11] = vec_i[2, :].dot(vec_j[3, :])
# angle16[12] = vec_i[3, :].dot(vec_j[0, :])
# angle16[13] = vec_i[3, :].dot(vec_j[1, :])
# angle16[14] = vec_i[3, :].dot(vec_j[2, :])
# angle16[15] = vec_i[3, :].dot(vec_j[3, :])
tmin = angle16.max() # arccos will be applied, so now swapped for max and min.
tmax = angle16.min()
adm[i, j] = tmin # i<j so minimum.
adm[j, i] = tmax
# Here adm is the inner product between unit vecrtos.
adm = np.where(adm > 1, 1, adm) # Numerically sometimes >1 value comes. In this case, convert to 1.
adm = np.where(adm < -1, -1, adm)
adm = np.arccos(adm) * 180. / np.pi
self.angular_distance_matrix = adm
def __str__(self):
s = "SimpleGridSphere nLon={}, nLat={}, nData={}.".format(self.nlon, self.nlat, self.get_counts().sum())
return s