apexpy.helpers.gc2gdlat() - Code Metrics - Inspection of "Merge pull request #59 from aburrell/rc_v1.1.0" - aburrell/apexpy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Push — main ( a2e919...ba9c59 )

by Angeline

created 2021-03-05 15:55 UTC

apexpy.helpers.gc2gdlat() A

↳ Parent: apexpy.helpers

Complexity

Conditions

Size

Total Lines	16
Code Lines	3

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	3
dl	0
loc	16
rs	10
c	0
b	0
f	0
cc	1
nop	1

# -*- coding: utf-8 -*-

"""This module contains helper functions used by :class:`~apexpy.Apex`."""

from __future__ import division, print_function, absolute_import

import time
import datetime as dt
import numpy as np


def checklat(lat, name='lat'):
    """Makes sure the latitude is inside [-90, 90], clipping close values
    (tolerance 1e-4).

    Parameters
    ----------
    lat : array-like
        latitude
    name : str, optional
        parameter name to use in the exception message

    Returns
    -------
    lat : ndarray or float
        Same as input where values just outside the range have been
        clipped to [-90, 90]

    Raises
    ------
    ValueError
        if any values are too far outside the range [-90, 90]
    """
    if np.any(np.abs(lat) > 90 + 1e-5):
        raise ValueError(name + ' must be in [-90, 90]')

    return np.clip(lat, -90.0, 90.0)


def getsinIm(alat):
    """Computes sinIm from modified apex latitude.

    Parameters
    ----------
    alat : array-like
        Modified apex latitude

    Returns
    -------
    sinIm : ndarray or float

    """

    alat = np.float64(alat)

    return 2 * np.sin(np.radians(alat)) / np.sqrt(4 - 3
                                                  * np.cos(np.radians(alat))**2)


def getcosIm(alat):
    """Computes cosIm from modified apex latitude.

    Parameters
    ----------
    alat : array-like
        Modified apex latitude

    Returns
    -------
    cosIm : ndarray or float

    """

    alat = np.float64(alat)

    return np.cos(np.radians(alat)) / np.sqrt(4 - 3
                                              * np.cos(np.radians(alat))**2)


def toYearFraction(date):
    """Converts :class:`datetime.date` or :class:`datetime.datetime` to decimal
    year.

    Parameters
    ----------
    date : :class:`datetime.date` or :class:`datetime.datetime`

    Returns
    -------
    year : float
        Decimal year

    Notes
    -----
    The algorithm is taken from http://stackoverflow.com/a/6451892/2978652

    """

    def sinceEpoch(date):
        """returns seconds since epoch"""
        return time.mktime(date.timetuple())
    year = date.year
    startOfThisYear = dt.datetime(year=year, month=1, day=1)
    startOfNextYear = dt.datetime(year=year + 1, month=1, day=1)

    yearElapsed = sinceEpoch(date) - sinceEpoch(startOfThisYear)
    yearDuration = sinceEpoch(startOfNextYear) - sinceEpoch(startOfThisYear)
    fraction = yearElapsed / yearDuration

    return date.year + fraction


def gc2gdlat(gclat):
    """Converts geocentric latitude to geodetic latitude using WGS84.

    Parameters
    ---------
    gclat : array-like
        Geocentric latitude

    Returns
    -------
    gdlat : ndarray or float
        Geodetic latitude

    """
    WGS84_e2 = 0.006694379990141317  # WGS84 first eccentricity squared
    return np.rad2deg(-np.arctan(np.tan(np.deg2rad(gclat)) / (WGS84_e2 - 1)))


def subsol(datetime):
    """Finds subsolar geocentric latitude and longitude.

    Parameters
    ----------
    datetime : :class:`datetime.datetime` or :class:`numpy.ndarray[datetime64]`
        Date and time in UTC (naive objects are treated as UTC)

    Returns
    -------
    sbsllat : float
        Latitude of subsolar point
    sbsllon : float
        Longitude of subsolar point

    Notes
    -----
    Based on formulas in Astronomical Almanac for the year 1996, p. C24.
    (U.S. Government Printing Office, 1994). Usable for years 1601-2100,
    inclusive. According to the Almanac, results are good to at least 0.01
    degree latitude and 0.025 degrees longitude between years 1950 and 2050.
    Accuracy for other years has not been tested. Every day is assumed to have
    exactly 86400 seconds; thus leap seconds that sometimes occur on December
    31 are ignored (their effect is below the accuracy threshold of the
    algorithm).

    After Fortran code by A. D. Richmond, NCAR. Translated from IDL
    by K. Laundal.

    """
    # Convert to year, day of year and seconds since midnight
    if isinstance(datetime, dt.datetime):
        year = np.asanyarray([datetime.year])
        doy = np.asanyarray([datetime.timetuple().tm_yday])
        ut = np.asanyarray([datetime.hour * 3600 + datetime.minute * 60
                            + datetime.second])
    elif isinstance(datetime, np.ndarray):
        # This conversion works for datetime of wrong precision or unit epoch
        times = datetime.astype('datetime64[s]')
        year_floor = times.astype('datetime64[Y]')
        day_floor = times.astype('datetime64[D]')
        year = year_floor.astype(int) + 1970
        doy = (day_floor - year_floor).astype(int) + 1
        ut = (times.astype('datetime64[s]') - day_floor).astype(float)
    else:
        raise ValueError("input must be datetime.datetime or numpy array")

    if not (np.all(1601 <= year) and np.all(year <= 2100)):
        raise ValueError('Year must be in [1601, 2100]')

    yr = year - 2000

    nleap = np.floor((year - 1601.0) / 4.0).astype(int)
    nleap -= 99
    mask_1900 = year <= 1900
    if np.any(mask_1900):
        ncent = np.floor((year[mask_1900] - 1601.0) / 100.0).astype(int)
        ncent = 3 - ncent
        nleap[mask_1900] = nleap[mask_1900] + ncent

    l0 = -79.549 + (-0.238699 * (yr - 4.0 * nleap) + 3.08514e-2 * nleap)
    g0 = -2.472 + (-0.2558905 * (yr - 4.0 * nleap) - 3.79617e-2 * nleap)

    # Days (including fraction) since 12 UT on January 1 of IYR:
    df = (ut / 86400.0 - 1.5) + doy

    # Mean longitude of Sun:
    lmean = l0 + 0.9856474 * df

    # Mean anomaly in radians:
    grad = np.radians(g0 + 0.9856003 * df)

    # Ecliptic longitude:
    lmrad = np.radians(lmean + 1.915 * np.sin(grad)
                       + 0.020 * np.sin(2.0 * grad))
    sinlm = np.sin(lmrad)

    # Obliquity of ecliptic in radians:
    epsrad = np.radians(23.439 - 4e-7 * (df + 365 * yr + nleap))

    # Right ascension:
    alpha = np.degrees(np.arctan2(np.cos(epsrad) * sinlm, np.cos(lmrad)))

    # Declination, which is also the subsolar latitude:
    sslat = np.degrees(np.arcsin(np.sin(epsrad) * sinlm))

    # Equation of time (degrees):
    etdeg = lmean - alpha
    nrot = np.round(etdeg / 360.0)
    etdeg = etdeg - 360.0 * nrot

    # Subsolar longitude calculation. Earth rotates one degree every 240 s.
    sslon = 180.0 - (ut / 240.0 + etdeg)
    nrot = np.round(sslon / 360.0)
    sslon = sslon - 360.0 * nrot

    # Return a single value from the output if the input was a single value
    if isinstance(datetime, dt.datetime):
        return sslat[0], sslon[0]
    return sslat, sslon


1			# -- coding: utf-8 --
2
3			"""This module contains helper functions used by :class:`~apexpy.Apex`."""
4
5			from __future__ import division, print_function, absolute_import
6
7			import time
8			import datetime as dt
9			import numpy as np
10
11
12			def checklat(lat, name='lat'):
13			"""Makes sure the latitude is inside [-90, 90], clipping close values
14			(tolerance 1e-4).
15
16			Parameters
17			----------
18			lat : array-like
19			latitude
20			name : str, optional
21			parameter name to use in the exception message
22
23			Returns
24			-------
25			lat : ndarray or float
26			Same as input where values just outside the range have been
27			clipped to [-90, 90]
28
29			Raises
30			------
31			ValueError
32			if any values are too far outside the range [-90, 90]
33			"""
34			if np.any(np.abs(lat) > 90 + 1e-5):
35			raise ValueError(name + ' must be in [-90, 90]')
36
37			return np.clip(lat, -90.0, 90.0)
38
39
40			def getsinIm(alat):
41			"""Computes sinIm from modified apex latitude.
42
43			Parameters
44			----------
45			alat : array-like
46			Modified apex latitude
47
48			Returns
49			-------
50			sinIm : ndarray or float
51
52			"""
53
54			alat = np.float64(alat)
55
56			return 2 * np.sin(np.radians(alat)) / np.sqrt(4 - 3
57			* np.cos(np.radians(alat))**2)
58
59
60			def getcosIm(alat):
61			"""Computes cosIm from modified apex latitude.
62
63			Parameters
64			----------
65			alat : array-like
66			Modified apex latitude
67
68			Returns
69			-------
70			cosIm : ndarray or float
71
72			"""
73
74			alat = np.float64(alat)
75
76			return np.cos(np.radians(alat)) / np.sqrt(4 - 3
77			* np.cos(np.radians(alat))**2)
78
79
80			def toYearFraction(date):
81			"""Converts :class:`datetime.date` or :class:`datetime.datetime` to decimal
82			year.
83
84			Parameters
85			----------
86			date : :class:`datetime.date` or :class:`datetime.datetime`
87
88			Returns
89			-------
90			year : float
91			Decimal year
92
93			Notes
94			-----
95			The algorithm is taken from http://stackoverflow.com/a/6451892/2978652
96
97			"""
98
99			def sinceEpoch(date):
100			"""returns seconds since epoch"""
101			return time.mktime(date.timetuple())
102			year = date.year
103			startOfThisYear = dt.datetime(year=year, month=1, day=1)
104			startOfNextYear = dt.datetime(year=year + 1, month=1, day=1)
105
106			yearElapsed = sinceEpoch(date) - sinceEpoch(startOfThisYear)
107			yearDuration = sinceEpoch(startOfNextYear) - sinceEpoch(startOfThisYear)
108			fraction = yearElapsed / yearDuration
109
110			return date.year + fraction
111
112
113			def gc2gdlat(gclat):
114			"""Converts geocentric latitude to geodetic latitude using WGS84.
115
116			Parameters
117			---------
118			gclat : array-like
119			Geocentric latitude
120
121			Returns
122			-------
123			gdlat : ndarray or float
124			Geodetic latitude
125
126			"""
127			WGS84_e2 = 0.006694379990141317 # WGS84 first eccentricity squared
128			return np.rad2deg(-np.arctan(np.tan(np.deg2rad(gclat)) / (WGS84_e2 - 1)))
129
130
131			def subsol(datetime):
132			"""Finds subsolar geocentric latitude and longitude.
133
134			Parameters
135			----------
136			datetime : :class:`datetime.datetime` or :class:`numpy.ndarray[datetime64]`
137			Date and time in UTC (naive objects are treated as UTC)
138
139			Returns
140			-------
141			sbsllat : float
142			Latitude of subsolar point
143			sbsllon : float
144			Longitude of subsolar point
145
146			Notes
147			-----
148			Based on formulas in Astronomical Almanac for the year 1996, p. C24.
149			(U.S. Government Printing Office, 1994). Usable for years 1601-2100,
150			inclusive. According to the Almanac, results are good to at least 0.01
151			degree latitude and 0.025 degrees longitude between years 1950 and 2050.
152			Accuracy for other years has not been tested. Every day is assumed to have
153			exactly 86400 seconds; thus leap seconds that sometimes occur on December
154			31 are ignored (their effect is below the accuracy threshold of the
155			algorithm).
156
157			After Fortran code by A. D. Richmond, NCAR. Translated from IDL
158			by K. Laundal.
159
160			"""
161			# Convert to year, day of year and seconds since midnight
162			if isinstance(datetime, dt.datetime):
163			year = np.asanyarray([datetime.year])
164			doy = np.asanyarray([datetime.timetuple().tm_yday])
165			ut = np.asanyarray([datetime.hour * 3600 + datetime.minute * 60
166			+ datetime.second])
167			elif isinstance(datetime, np.ndarray):
168			# This conversion works for datetime of wrong precision or unit epoch
169			times = datetime.astype('datetime64[s]')
170			year_floor = times.astype('datetime64[Y]')
171			day_floor = times.astype('datetime64[D]')
172			year = year_floor.astype(int) + 1970
173			doy = (day_floor - year_floor).astype(int) + 1
174			ut = (times.astype('datetime64[s]') - day_floor).astype(float)
175			else:
176			raise ValueError("input must be datetime.datetime or numpy array")
177
178			if not (np.all(1601 <= year) and np.all(year <= 2100)):
179			raise ValueError('Year must be in [1601, 2100]')
180
181			yr = year - 2000
182
183			nleap = np.floor((year - 1601.0) / 4.0).astype(int)
184			nleap -= 99
185			mask_1900 = year <= 1900
186			if np.any(mask_1900):
187			ncent = np.floor((year[mask_1900] - 1601.0) / 100.0).astype(int)
188			ncent = 3 - ncent
189			nleap[mask_1900] = nleap[mask_1900] + ncent
190
191			l0 = -79.549 + (-0.238699 * (yr - 4.0 * nleap) + 3.08514e-2 * nleap)
192			g0 = -2.472 + (-0.2558905 * (yr - 4.0 * nleap) - 3.79617e-2 * nleap)
193
194			# Days (including fraction) since 12 UT on January 1 of IYR:
195			df = (ut / 86400.0 - 1.5) + doy
196
197			# Mean longitude of Sun:
198			lmean = l0 + 0.9856474 * df
199
200			# Mean anomaly in radians:
201			grad = np.radians(g0 + 0.9856003 * df)
202
203			# Ecliptic longitude:
204			lmrad = np.radians(lmean + 1.915 * np.sin(grad)
205			+ 0.020 * np.sin(2.0 * grad))
206			sinlm = np.sin(lmrad)
207
208			# Obliquity of ecliptic in radians:
209			epsrad = np.radians(23.439 - 4e-7 * (df + 365 * yr + nleap))
210
211			# Right ascension:
212			alpha = np.degrees(np.arctan2(np.cos(epsrad) * sinlm, np.cos(lmrad)))
213
214			# Declination, which is also the subsolar latitude:
215			sslat = np.degrees(np.arcsin(np.sin(epsrad) * sinlm))
216
217			# Equation of time (degrees):
218			etdeg = lmean - alpha
219			nrot = np.round(etdeg / 360.0)
220			etdeg = etdeg - 360.0 * nrot
221
222			# Subsolar longitude calculation. Earth rotates one degree every 240 s.
223			sslon = 180.0 - (ut / 240.0 + etdeg)
224			nrot = np.round(sslon / 360.0)
225			sslon = sslon - 360.0 * nrot
226
227			# Return a single value from the output if the input was a single value
228			if isinstance(datetime, dt.datetime):
229			return sslat[0], sslon[0]
230			return sslat, sslon
231

aburrell / apexpy

Push — main ( a2e919...ba9c59 )

apexpy.helpers.gc2gdlat() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like