The package can be installed from r-universe by:
install.packages("LakeGeoid", repos = c("https://cavios.r-universe.dev", "https://cloud.r-project.org"))
Please find a description of all functions here
The package contain 3 test data sets
library(LakeGeoid)
library(sf)
data(kivu) #SWOT Raster data (250m) for lake Kivu
data(myshape) # shapefile for lake kivu
data(C2kivu) # CryoSat-2 data for lake Kivu
#myshape <- st_read("myshape.shp")
#simplify shapefile to reduce the amount of triangles in the mesh
shpsimple<-st_simplify(myshape,dTolerance=300)
fit<-getLakeGeoid(kivu,maxEdge=3000,myshape=shpsimple,UTM=TRUE)
summary(fit)
-------------------------
Summary for getLakeGeoid
-------------------------
Converged with a negative log likelihood of -177695.689
Number of parameters: 4
Par 1 : sdObs = 0.072
Par 2 : Tau = 1100142845.643
Par 3 : Kappa = 0
Par 4 : SigmaRW = 2.639
Plotting the model output
plot(fit)
We can get geoid correction at the coordinates of CryoSat-2
C2new<-getGeoidCorr(fit,C2kivu)
head(C2new)
time lat lon height satid ocval sat gCorr
1 2010.559 -2.332386 29.12126 1462.901 9 97 C2E 0.1911031
2 2010.559 -2.329602 29.12361 1462.984 9 97 C2E 0.1859167
3 2010.559 -2.326954 29.12208 1462.767 9 98 C2E 0.1833408
4 2010.559 -2.318522 29.13137 1462.518 9 98 C2E 0.1705378
5 2010.559 -2.315892 29.12931 1461.234 9 97 C2E 0.1637882
6 2010.559 -2.313156 29.13028 1463.349 9 97 C2E 0.1558789
The column gCorr contains the geoid correction which sould be added to the heights
SWOT Level 2 Water Mask Raster Image 250m data can be downloaded directly from PO.DAAC for a lake
in the SWOT Prior Lake Database (PLD). This requires a (free) NASA Earthdata account.
Add your login to the file ~/.netrc:
machine urs.earthdata.nasa.gov login <username> password <password>
(alternatively set the environment variables EARTHDATA_USER and EARTHDATA_PASS).
library(LakeGeoid)
library(sf)
pld <- "SWOTPLDv2.2_simple/lakesub_simple.shp" # PLD shapefile
lake <- getPLDLake(1320035532, pld) # lake Kivu
dat <- getSWOTRaster(lake, startDate = "2024-01-01", endDate = "2024-06-30",
maxDarkFrac = 0.1, outDir = "SWOT", keep = TRUE)
fit <- getLakeGeoid(dat, maxEdge = 3000, myshape = st_simplify(lake, dTolerance = 300))
The lake can also be given as a polygon of the class "sf" (e.g. getSWOTRaster(myshape)) or as a lake id
together with the PLD (getSWOTRaster(1320035532, pld = pld)). By default version D of the product is used;
use shortName = "SWOT_L2_HR_Raster_250m_2.0" for version C.
flagOutliers() adds a column flag with "ok" or the reason an observation is flagged:
dark water (darkwater_frac), optional thresholds on wse_uncert, N_PIXC and distance to shore,
overpasses with too few observations ("sparse") or a large tilt relative to the mean field ("tilt"),
and remaining gross errors relative to the overpass median and the mean field ("mad").
A table with statistics per overpass is attached as attr(dat, "overpass").
dat <- flagOutliers(dat)
plotOverpass(dat, lake, cycle = 5, pass = 249) # flagged observations shown as gray crosses
fit <- getLakeGeoid(dat[dat$flag == "ok", ], maxEdge = 3000,
myshape = st_simplify(lake, dTolerance = 300), cellSize = 750)
SWOT passes can have small constant offsets. With passBias = TRUE a bias per pass is estimated
(the biases sum to zero); the estimates are in fit$passBias and shown by summary(fit).
The estimated field can be exported as a raster (GeoTIFF) on a regular grid:
fit <- getLakeGeoid(dat[dat$flag == "ok", ], maxEdge = 3000, myshape = shp, cellSize = 750, passBias = TRUE)
# optionally also a random east/north tilt per overpass (e.g. SWOT roll errors); estimates in fit$tilt
fit <- getLakeGeoid(dat[dat$flag == "ok", ], maxEdge = 3000, myshape = shp, cellSize = 750,
passBias = TRUE, overpassTilt = TRUE)
plot(fit) # use plot(fit, showMesh = TRUE) to see the triangles
r <- getGeoidGrid(fit, res = 250)
terra::writeRaster(r, "geoidcorr.tif")