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

Completed

Push — main ( fcb71d...fb1ae8 )

by Angeline

created 2025-01-07 15:08 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`."""

import datetime as dt
import numpy as np
import time


def set_array_float(in_val):
    """Set array data type to float.

    Parameters
    ----------
    in_val : any
        Input value, only modified if it is a np.ndarray

    Returns
    -------
    out_val : any
        Output value, if `in_val` was an array, `out_val` will be an array of
        type `np.float64`.

    """

    if isinstance(in_val, np.ndarray):
        out_val = in_val.astype(np.float64)
    else:
        out_val = in_val

    return out_val


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('{:} must be in [-90, 90]'.format(name))

    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`
        Input date or datetime object

    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

    year += fraction

    return year


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

aburrell / apexpy

Push — main ( fcb71d...fb1ae8 )

apexpy.helpers.gc2gdlat() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like