harphub / meteogrid

R package for working with gridded meteorological data
https://harphub.github.io/meteogrid/
Other
0 stars 6 forks source link

Bug: Incorrect mapfactor calculation for .geowind.LCC #21

Open kobebryant432 opened 3 weeks ago

kobebryant432 commented 3 weeks ago

There seems to be an error in the calculation of the mapfactor for LCC projection

https://github.com/harphub/meteogrid/blob/9da2345441037e82a518d4e7a18145285ae9a5e2/R/geowind.R#L177C1-L178C107

In particular the power of the first term should be $$1-sin(\phi_1)$$ as apposed to $$1-cos(\phi_1)$$

  mapfactor <- (refcos/cos(lalo$lat * rad))^(1 - refsin) * ((1 + refsin)/(1 + sin(lalo$lat * rad)))^refsin
  angle <- refsin * (lalo$lon - reflon) * rad

Credits for finding this -Wout Dewettinck

@adeckmyn

kobebryant432 commented 3 weeks ago

Derivation

Formulas

We define the coordinates in the LCC $xy$-system as follows: $x(\lambda, \phi) = \rho \sin \left(\sin \phi_1 (\lambda - \lambda_0)  \right)$ $y(\lambda, \phi) = \rho_0 - \rho \cos \left(\sin \phi_1 (\lambda - \lambda_0)  \right)$ with $\phi_1$ the reference latitude (0° at the equator, 90° at the North Pole) and

$\rho(\phi) = a \frac{\cos \phi_1}{\sin \phi_1} \left(\frac{\tan \left(\frac{\pi}{4} + \frac{\phi_1}{2} \right)}{\tan \left(\frac{\pi}{4} + \frac{\phi}{2} \right)} \right)^{\sin \phi_1}$

The value $\rho_0$ is defined such that $y = 0$ at the equator ($\phi$ = 0):

$\rho_0 = a \frac{\cos \phi_1}{\sin \phi_1} \tan \left(\frac{\pi}{4} + \frac{\phi_1}{2} \right)^{\sin \phi_1}$

The value $a$ is the radius of the Earth.

Calculating the map factor (Cylindrical projection) $h$ - scale factor along meridian of longitude $k$ - scale factor along parallel of latitude

For the LCC this is equivalent to from https://pubs.usgs.gov/pp/1395/report.pdf formula (15-4) $h = k = \frac{\rho}{R\cos{\phi}} = \cos \phi_1 \tan^{\sin \phi_1} \left( \frac{\pi}{4} + \frac{\phi_1}{2} \right) \Big/ \left[ \cos \phi \tan^{\sin \phi_1} \left( \frac{\pi}{4} + \frac{\phi}{2} \right) \right]$

using the fact that $\tan \left( \frac{\pi}{4} + \frac{x}{2} \right) = \frac{1 + \sin x}{\cos x}$ gives:

$k=\left( \frac{\cos \phi_1 }{\cos \phi } \right)^{1 - \sin \phi_1 } \times \left( \frac{1 + \sin \phi_1 }{1 + \sin \phi } \right)^{\sin \phi_1 }$