Skip to content

Commit

Permalink
Improve harmonic tide calculations
Browse files Browse the repository at this point in the history
  • Loading branch information
@kylecorry31
kylecorry31 committed Aug 1, 2024
1 parent d98c678 commit c9b3692
Show file tree
Hide file tree
Showing 6 changed files with 155 additions and 84 deletions.
2 changes: 1 addition & 1 deletion build.gradle.kts
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ plugins {
}

group = "com.kylecorry"
version = "9.8.0"
version = "9.8.1"

afterEvaluate {
publishing {
Expand Down
46 changes: 23 additions & 23 deletions src/main/kotlin/com/kylecorry/sol/science/oceanography/AstronomicalArgumentCalculator.kt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,27 +6,7 @@ internal object AstronomicalArgumentCalculator {

fun get(constituent: TideConstituent, year: Int): Float {
// From https://www.dfo-mpo.gc.ca/science/data-donnees/tidal-marees/argument-u-v-eng.html
return if (year <= 2021) {
get2021(constituent)
} else {
get2022(constituent)
}
}

private fun get2021(constituent: TideConstituent): Float {
return when (constituent) {
TideConstituent.M2 -> 304.098f
TideConstituent.S2 -> 0.138f
TideConstituent.N2 -> 33.743f
TideConstituent.K1 -> 2.518f
TideConstituent.M4 -> 248.196f
TideConstituent.O1 -> 304.678f
TideConstituent.P1 -> 348.423f
TideConstituent.L2 -> 17.969f // Calculated
TideConstituent.K2 -> 184.784f
TideConstituent.MS4 -> 304.236f
TideConstituent.Z0 -> 0f
}
return get2022(constituent)
}

private fun get2022(constituent: TideConstituent): Float {
Expand All @@ -35,12 +15,32 @@ internal object AstronomicalArgumentCalculator {
TideConstituent.S2 -> 0.113f
TideConstituent.N2 -> 46.36f
TideConstituent.K1 -> 3.577f
TideConstituent.M4 -> 90.026f // M2 * 2
TideConstituent.M4 -> 45.013f * 2
TideConstituent.O1 -> 44.031f
TideConstituent.M6 -> 3 * 45.013f
TideConstituent.MK3 -> 45.013f + 3.577f // M2 + K1
TideConstituent.S4 -> 0.113f * 2 // S2 + S2
TideConstituent.MN4 -> 45.013f + 46.36f // M2 + N2
TideConstituent.NU2 -> 45.013f // M2
TideConstituent.S6 -> 0.113f * 3 // S2 + S2 + S2
TideConstituent.MU2 -> 45.013f * 2 - 0.113f // M2 * 2 - S2
TideConstituent._2N2 -> 45.013f // M2
TideConstituent.LAM2 -> 45.013f // M2
TideConstituent.S1 -> 0f
TideConstituent.SSA -> 0f
TideConstituent.SA -> 0f
TideConstituent.MSF -> 0.113f - 45.013f // S2 - M2
TideConstituent.RHO -> 45.013f - 3.577f // NU2 - K1
TideConstituent.Q1 -> 45.003f
TideConstituent.T2 -> 0f
TideConstituent.R2 -> 0f
TideConstituent.P1 -> 348.805f
TideConstituent._2SM2 -> 2 * 0.113f - 45.013f // 2 * S2 - M2
TideConstituent.L2 -> 213.7788f // Calculated
TideConstituent._2MK3 -> 45.013f + 44.031f // M2 + O1
TideConstituent.K2 -> 186.65f
TideConstituent.MS4 -> 45.126f // M2 + S2
TideConstituent.M8 -> 45.013f * 4 // M2 * 4
TideConstituent.MS4 -> 45.013f + 0.113f // M2 + S2
TideConstituent.Z0 -> 0f
}
}
Expand Down
44 changes: 23 additions & 21 deletions src/main/kotlin/com/kylecorry/sol/science/oceanography/NodeFactorCalculator.kt
View file
Original file line number Diff line number Diff line change
@@ -1,33 +1,15 @@
package com.kylecorry.sol.science.oceanography

import com.kylecorry.sol.math.SolMath.cube
import com.kylecorry.sol.math.SolMath.power
import com.kylecorry.sol.math.SolMath.square

internal object NodeFactorCalculator {
// TODO: Use Schureman equations to calculate (https://github.com/sam-cox/pytides/blob/master/pytides/nodal_corrections.py)

fun get(constituent: TideConstituent, year: Int): Float {
// From https://www.dfo-mpo.gc.ca/science/data-donnees/tidal-marees/facteur-node-factor-eng.html
return if (year <= 2021) {
get2021(constituent)
} else {
get2022(constituent)
}
}

private fun get2021(constituent: TideConstituent): Float {
return when (constituent) {
TideConstituent.M2 -> 0.986f
TideConstituent.S2 -> 1.001f
TideConstituent.N2 -> 0.984f
TideConstituent.K1 -> 1.053f
TideConstituent.M4 -> square(0.986f)
TideConstituent.O1 -> 1.088f
TideConstituent.P1 -> 0.997f
TideConstituent.L2 -> 0.8537f // Calculated
TideConstituent.K2 -> 1.125f
TideConstituent.MS4 -> 0.986f * 1.001f // M2 * S2
TideConstituent.Z0 -> 1f
}
return get2022(constituent)
}

private fun get2022(constituent: TideConstituent): Float {
Expand All @@ -38,9 +20,29 @@ internal object NodeFactorCalculator {
TideConstituent.K1 -> 1.081f
TideConstituent.M4 -> square(0.975f)
TideConstituent.O1 -> 1.140f
TideConstituent.M6 -> cube(0.975f)
TideConstituent.MK3 -> 0.975f * 1.081f // M2 * K1
TideConstituent.S4 -> square(1.001f) // S2 * S2
TideConstituent.MN4 -> 0.975f * 0.978f // M2 * N2
TideConstituent.NU2 -> 0.975f // M2
TideConstituent.S6 -> cube(1.001f) // S2 * S2 * S2
TideConstituent.MU2 -> square(0.975f) * 1.001f // M2 * M2 * S2
TideConstituent._2N2 -> 0.975f // M2
TideConstituent.LAM2 -> 0.975f // M2
TideConstituent.S1 -> 1f
TideConstituent.SSA -> 1f
TideConstituent.SA -> 1f
TideConstituent.MSF -> 1.001f * 0.975f // S2 * M2
TideConstituent.RHO -> 0.975f * 1.081f // NU2 * K1
TideConstituent.Q1 -> 1.102f
TideConstituent.T2 -> 1f
TideConstituent.R2 -> 1f
TideConstituent.P1 -> 0.995f
TideConstituent._2SM2 -> square(1.001f) * 0.975f // 2 * S2 * M2
TideConstituent.L2 -> 1.2437f
TideConstituent._2MK3 -> 0.975f * 1.140f // M2 * O1
TideConstituent.K2 -> 1.214f
TideConstituent.M8 -> power(0.975f, 4) // M2 * M2 * M2 * M2
TideConstituent.MS4 -> 0.975f * 1.001f
TideConstituent.Z0 -> 1f
}
Expand Down
27 changes: 27 additions & 0 deletions src/main/kotlin/com/kylecorry/sol/science/oceanography/TidalHarmonic.kt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,36 @@ enum class TideConstituent(val id: Long, val speed: Float) {
K1(4, 15.041069f),
M4(5, 57.96821f),
O1(6, 13.943035f),
M6(7, 86.95232f),
MK3(8, 44.025173f),
S4(9, 60f),
MN4(10, 57.423832f),
NU2(11, 28.512583f),
S6(12, 90f),
MU2(13, 27.968208f),
_2N2(14, 27.895355f),
// OO1(15, 16.139101f),
LAM2(16, 29.455626f),
S1(17, 15f),
// M1(18, 14.496694f),
// J1(19, 15.5854435f),
// MM(20, 0.5443747f),
SSA(21, 0.0821373f),
SA(22, 0.0410686f),
MSF(23, 1.0158958f),
// MF(24, 1.0980331f),
RHO(25, 13.471515f),
Q1(26, 13.398661f),
T2(27, 29.958933f),
R2(28, 30.041067f),
// _2Q1(29, 12.854286f),
P1(30, 14.958931f),
_2SM2(31, 31.015896f),
// M3(32, 43.47616f),
L2(33, 29.528479f),
_2MK3(34, 42.92714f),
K2(35, 30.082138f),
M8(36, 115.93642f),
MS4(37, 58.984104f),
Z0(0, 0f)
}
6 changes: 1 addition & 5 deletions .../kotlin/com/kylecorry/sol/science/oceanography/waterlevel/HarmonicWaterLevelCalculator.kt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,11 +15,7 @@ class HarmonicWaterLevelCalculator(private val harmonics: List<TidalHarmonic>) :
IWaterLevelCalculator {
override fun calculate(time: ZonedDateTime): Float {
// TODO: This is a hack to get the correct year for the astronomical arguments
val year = if (time.year <= 2021){
2021
} else {
2022
}
val year = 2022
val start = ZonedDateTime.of(
LocalDateTime.of(year, 1, 1, 0, 0),
ZoneId.of("UTC")
Expand Down
114 changes: 80 additions & 34 deletions src/test/kotlin/com/kylecorry/sol/science/oceanography/OceanographyServiceTest.kt
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,14 @@
package com.kylecorry.sol.science.oceanography

import com.kylecorry.sol.math.SolMath
import com.kylecorry.sol.math.SolMath.roundPlaces
import com.kylecorry.sol.math.SolMath.toDegrees
import com.kylecorry.sol.math.statistics.Statistics
import com.kylecorry.sol.science.astronomy.locators.Moon
import com.kylecorry.sol.science.oceanography.waterlevel.*
import com.kylecorry.sol.time.Time
import com.kylecorry.sol.time.Time.atEndOfDay
import com.kylecorry.sol.time.Time.atStartOfDay
import com.kylecorry.sol.time.Time.roundNearestMinute
import com.kylecorry.sol.units.*
import org.junit.jupiter.api.Assertions.assertEquals
Expand All @@ -15,6 +19,7 @@ import java.io.File
import java.net.http.HttpClient
import java.time.*
import java.util.stream.Stream
import kotlin.math.acos

internal class OceanographyServiceTest {
@Test
Expand All @@ -37,44 +42,85 @@ internal class OceanographyServiceTest {

@Test
fun getTidesHarmonics() {
val harmonics = listOf(
TidalHarmonic(TideConstituent.M2, 1.66f, 2.3f),
TidalHarmonic(TideConstituent.S2, 0.35f, 25f),
TidalHarmonic(TideConstituent.N2, 0.41f, 345.8f),
TidalHarmonic(TideConstituent.K1, 0.2f, 166.1f),
TidalHarmonic(TideConstituent.O1, 0.15f, 202f),
TidalHarmonic(TideConstituent.P1, 0.07f, 176.6f),
TidalHarmonic(TideConstituent.M4, 0.19f, 35.8f),
TidalHarmonic(TideConstituent.K2, 0.1f, 21.7f),
TidalHarmonic(TideConstituent.L2, 0.04f, 349.9f),
TidalHarmonic(TideConstituent.MS4, 0.05f, 106.4f),
TidalHarmonic(TideConstituent.Z0, 0f, 0f)
)

val zone = ZoneId.of("America/New_York")
val date = LocalDate.of(2023, 12, 22).atStartOfDay(zone)

val expected = listOf(
Tide.high(ZonedDateTime.of(date.toLocalDate(), LocalTime.of(3, 31), zone)),
Tide.low(ZonedDateTime.of(date.toLocalDate(), LocalTime.of(10, 23), zone)),
Tide.high(ZonedDateTime.of(date.toLocalDate(), LocalTime.of(15, 56), zone)),
Tide.low(ZonedDateTime.of(date.toLocalDate(), LocalTime.of(21, 36), zone))
)

val service = OceanographyService()
val constituents = listOf(
Triple(1, 0.505968f, 2.3f),
Triple(2, 0.10668f, 25.0f),
Triple(3, 0.124968f, 345.8f),
Triple(4, 0.06096000000000001f, 166.1f),
Triple(5, 0.057912000000000005f, 35.8f),
Triple(6, 0.045720000000000004f, 202.0f),
Triple(7, 0.006096000000000001f, 220.1f),
Triple(8, 0.009144f, 19.5f),
Triple(9, 0.006096000000000001f, 5.1f),
Triple(10, 0.027432f, 347.9f),
Triple(11, 0.021336000000000004f, 341.0f),
Triple(12, 0.0f, 222.9f),
Triple(13, 0.024384000000000003f, 344.5f),
Triple(14, 0.018288f, 333.0f),
Triple(15, 0.0030480000000000004f, 195.7f),
Triple(16, 0.0030480000000000004f, 345.5f),
Triple(17, 0.006096000000000001f, 121.9f),
Triple(18, 0.0030480000000000004f, 204.0f),
Triple(19, 0.006096000000000001f, 181.1f),
Triple(20, 0.018288f, 73.9f),
Triple(21, 0.015240000000000002f, 75.1f),
Triple(22, 0.06096000000000001f, 145.3f),
Triple(23, 0.0f, 0.0f),
Triple(24, 0.0f, 0.0f),
Triple(25, 0.0030480000000000004f, 230.9f),
Triple(26, 0.012192000000000001f, 185.3f),
Triple(27, 0.009144f, 9.0f),
Triple(28, 0.0f, 252.2f),
Triple(29, 0.0030480000000000004f, 172.9f),
Triple(30, 0.021336000000000004f, 176.6f),
Triple(31, 0.0030480000000000004f, 48.8f),
Triple(32, 0.006096000000000001f, 34.1f),
Triple(33, 0.012192000000000001f, 349.9f),
Triple(34, 0.009144f, 357.1f),
Triple(35, 0.030480000000000004f, 21.7f),
Triple(36, 0.0f, 330.4f),
Triple(37, 0.015240000000000002f, 106.4f)
).mapNotNull {
val (order, amplitude, phase) = it
val constituent = TideConstituent.entries.firstOrNull { it.id == order.toLong() }
if (constituent == null) {
null
} else {
TidalHarmonic(constituent, amplitude, phase)
}
}

val tides = service.getTides(
HarmonicWaterLevelCalculator(harmonics),
date,
date.atEndOfDay()
)
val start = ZonedDateTime.of(2024, 7, 1, 0, 0, 0, 0, ZoneId.of("America/New_York"))
val end = ZonedDateTime.of(2024, 7, 29, 0, 0, 0, 0, ZoneId.of("America/New_York"))
val station = "8452660"
val url =
"https://api.tidesandcurrents.noaa.gov/api/prod/datagetter?product=predictions&application=NOS.COOPS.TAC.WL&begin_date=20240701&end_date=20240730&datum=MLLW&station=$station&time_zone=lst_ldt&units=english&interval=hilo&format=csv"
File("temp").mkdirs()

assertEquals(expected.size, tides.size)
val sourceString = if (File("temp/${station}.csv").exists()) {
File("temp/${station}.csv").readText()
} else {
// Download the file
val httpClient = HttpClient.newHttpClient()
val request = java.net.http.HttpRequest.newBuilder()
.uri(java.net.URI.create(url))
.build()
val response = httpClient.send(request, java.net.http.HttpResponse.BodyHandlers.ofString())
val body = response.body()
File("temp/${station}.csv").writeText(body)
body
}

for (i in tides.indices) {
assertEquals(expected[i].isHigh, tides[i].isHigh)
timeEquals(tides[i].time, expected[i].time, Duration.ofMinutes(30))
val source = sourceString.trim().lines().drop(1).map {
val parts = it.split(",")
Tide(
LocalDateTime.parse(parts[0].replace(" ", "T")).atZone(ZoneId.of("America/New_York")),
parts[2] == "H",
parts[1].toFloat()
)
}

testCalculator("Harmonic", HarmonicWaterLevelCalculator(constituents), start, end, source)
}

@Test
Expand Down

0 comments on commit c9b3692

Please sign in to comment.