About This Document

This document is an Emacs Org Mode plain-text file with code and text embedded. If you are viewing:

A DOC or PDF file, then it was generated by exporting from Org. Not all of the Org parts (code, results, comments, etc.) were exported. The Org source file is available upon request, and may be embedded in the PDF. Most non-Apple PDF viewers provide easy access to embedded or attached files.
A file with a org extension in something other than Emacs, then you are seeing the canonical version and the full source, but without any syntax highlighting, document structure, or the ability to execute the code blocks.
An Org file within Emacs, then this is the canonical version. You should be able to fully interact and reproduce the contents of this document, although it may require 3rd-party applications (Python, etc.) and a similar Emacs configuration. This is available upon request.

Workflow

To recreate this work

See the hacking.org file

After updates, re-run make, and then…

Run the workflow-update block below
- Cleaning all result blocks with C-u C-c C-v k or (org-babel-remove-result-one-or-many), then
- Executing all blocks (without :eval no) using C-c C-v C-b or (org-babel-execute-buffer)
Review and commit changes
Re-run the workflow-update so that exported files have the right git commit
- Review changes - there should be NONE
Push updates
- git push
- Upload data to dataverse

Summary

We have produced an open and reproducible estimate of Greenland ice sheet solid ice discharge from 1986 through last month. Our results show three modes at the total ice-sheet scale: Steady discharge from 1986 through 2000, increasing discharge from 2000 through 2005, steady discharge from 2005 through the recent record. The behavior of individual sectors and glaciers is more complicated. This work was done to provide a 100% reproducible and operational estimate to help constrain mass balance and sea level rise estimates.

New in this version

PROMICE ice velocity from Sentinel1 citep:solgaard_2021,solgaard_2021_data is now generated at 200 m rather than 500 m resolution. There is no significant code change associated with this updated data. This results in higher estimated velocities from the smallest glaciers, with a slight increase in overall discharge, and a more than slight increase in the SE sector due to the many small outlet glaciers there.

BedMachine has been updated to v5. Uncertainty has declined by ~10 % (that is, from ~10 % of total discharge to ~9 %, or from ~50 Gt yr^-1 to ~45 Gt yr^-1). With the BedMachine v5 update, the southeast sector is no longer patched with citet:millan_2018 data.

The elevation change code and data have been completely rewritten. Previous versions used the GIMP DEM (NSIDC 0645; citet:howat_2014_greenland,NSIDC_0645), a DEM generated from several years of data but with each pixel time-stamped. Elevation change was generated per-pixel using the citet:khan_2016_geodetic elevation change data set. The new workflow uses the PROMICE PRODEM 2020 DEM (CITE winstrup) as the baseline DEM, and annual PRODEMs forward in time from there. When no PRODEM exists (at the moment, for the past year) there is no computed change in elevation, and backward one year to 2019 (the first PRODEM year). Prior to 2019, we still use citet:khan_2016_geodetic to the beginning of that dataset which is 1995. Prior to 1995 there is no computed change in elevation. Annual DEMs and change in elevations are interpolated linearly and used to estimate a thickness at each velocity timestamp.

The majority of other velocity products from the NSIDC have also been updated. Current versions are captured in the metadata and digital workbooks published with this document. There are no significant changes with these updates which generally include recent data and occasionally reprocess historical data.

A new subsection titled “Limitations” within the Discussion section addresses some limitations of this product and warns against specific use-cases.

Introduction

The mass of the Greenland Ice Sheet is decreasing (e.g. citet:fettweis_2017_reconstructions,van-den-broeke_2017_greenland,wiese_2016_jpl,khan_2016_geodetic). Most ice sheet mass loss – as iceberg discharge, submarine melting, and meltwater runoff – enters the fjords and coastal seas, and therefore ice sheet mass loss directly contributes to sea-level rise citep:wcrp_2018,moon_2018_rising,nerem_2018_climate,chen_2017_increasing. Greenland’s total ice loss can be estimated through a variety of independent methods, for example direct mass change estimates from GRACE citep:wiese_2016_jpl or by using satellite altimetry to estimate surface elevation change, which is then converted into mass change (using a firn model, e.g., citet:khan_2016_geodetic). However, partitioning the mass loss between ice discharge (D) and surface mass balance (SMB) remains challenging (cf. citet:rignot_2008_mass and citet:enderlin_2014_improved). Correctly assessing mass loss, as well as the attribution of this loss (SMB or D), is critical to understanding the process-level response of the Greenland Ice Sheet to climate change and thus improving models of future ice sheet changes and associated sea-level rise citep:moon_2018_rising.

The total mass of an ice sheet, or a drainage basin, changes if the mass gain (SMB inputs, primarily snowfall) is not balanced by the mass loss (D and SMB outputs, the latter generally meltwater runoff). This change is typically termed ice sheet mass balance (MB) and the formal expression for this rate of change in mass is (e.g., citet:cuffey_2010_the-physics),

\begin{equation} \frac{\mathrm{d}M}{\mathrm{d}t} = ρ ∫_A b \, \mathrm{d}A - ∫_g Q \, \mathrm{d}g, \end{equation}

where $ρ$ is the average density of ice, $b$ is an area mass balance, and $Q$ is the discharge flux. The left-hand side of the equation is the rate of change of mass, the first term on the right-hand side is the area $A$ integrated surface mass balance (SMB), and the second term is the discharge $D$ mass flow rate that drains through gate $g$. Equation eq:dMdt is often simplified to

\begin{equation} MB = SMB - D \end{equation}

where $MB$ is the mass balance, and referred to as the “input–output” method (e.g., citet:khan_2015_greenland). Virtually all studies agree on the trend of Greenland mass balance, but large discrepancies persist in both the magnitude and attribution. Magnitude discrepancies include, for example, citet:kjeldsen_2015_spatial reporting a mass imbalance of -250 $±$ 21 Gt yr^-1 during 2003 to 2010, citet:ewert_2012_volume reporting -181 $±$ 28 Gt yr^-1 during 2003 to 2008, and citet:rignot_2008_mass reporting a mass imbalance of -265 $±$ 19 Gt yr^-1 during 2004 to 2008. Some of these differences may be due to different ice sheet area masks used in the studies. Attribution discrepancies include, for example, citet:enderlin_2014_improved attributing the majority (64 %) of mass loss to changes in SMB during the 2005 to 2009 period but citet:rignot_2008_mass attributing the majority (85 %) of mass loss to changes in D during the 2004 to 2008 period.

Discharge may be calculated through several methods, including mass flow rate through gates (e.g. citet:enderlin_2014_improved,king_2018_seasonal,mouginot_2019_forty), or solving as a residual from independent mass balance terms (e.g. citet:kjaer_2012_aerial,kjeldsen_2015_spatial). The gate method that we use in this study incorporates ice thickness and an estimated vertical profile from the observed surface velocity to calculate the discharge. A typical formulation of discharge across a gate $D_g$ is,

\begin{equation} D_g = ρ \, V \, H \, w, \end{equation}

where $ρ$ is the average density of ice, $V$ is depth-average gate-perpendicular velocity, $H$ is the ice thickness, and $w$ is the gate width. Uncertainties in $V$ and $H$ naturally influence the estimated discharge. At fast-flowing outlet glaciers, $V$ is typically assumed to be equal at all ice depths, and observed surface velocities can be directly translated into depth-averaged velocities (as in citet:enderlin_2014_improved,king_2018_seasonal). To minimize uncertainty from SMB or basal mass balance corrections downstream of a flux gate, the gate should be at the grounding line of the outlet glacier. Unfortunately, uncertainty in bed elevation (translating to ice thickness uncertainty) increases toward the grounding line.

Conventional methods of gate selection involve handpicking gate locations, generally as linear features (e.g., citet:enderlin_2014_improved) or visually approximating ice-orthogonal gates at one point in time (e.g., citet:king_2018_seasonal). Manual gate definition is suboptimal. For example, the largest discharging glaciers draw from an upstream radially diffusing region that may not easily be represented by a single linear gate. Approximately flow-orthogonal curved gates may not be flow orthogonal on the multidecade timescale due to changing flow directions. Manual gate selection makes it difficult to update gate locations, corresponding with glacier termini retreat or advance, in a systematic and reproducible fashion. We therefore adopt an algorithmic approach to generate gates based on a range of criteria.

Here, we present a discharge dataset based on gates selected in a reproducible fashion by a new algorithm. Relative to previous studies, we employ ice velocity observation over a longer period with higher temporal frequency and denser spatial coverage. We use ice velocity from 1986 through last month including twelve-day velocities for the last ~500 days of the time series, and discharge at 200 m pixel resolution capturing all ice flowing faster than 100 m yr^-1 that crosses glacier termini into fjords.

Input data

Historically, discharge gates were selected along well-constrained flight lines of airborne radar data citep:enderlin_2014_improved. Recent advances in ice thickness estimates through NASA Operation IceBridge citep:millan_2018_vulnerability, NASA Oceans Melting Greenland (OMG; citet:fenty_2016_oceans), fjord bathymetry citep:tinto_2015_bathymetry, and methods to estimate thickness from surface properties (e.g., citet:mcnabb_2012_using,james_2016_automated) have been combined into digital bed elevation models such as BedMachine citep:morlighem_2017_bedmachine,NSIDC_BedMachine or released as independent datasets citep:millan_2018_vulnerability. From these advances, digital bed elevation models have become more robust at tidewater glacier termini and grounding lines. The incorporation of flight-line ice thickness data into higher-level products that include additional methods and data means gates are no longer limited to flight lines (e.g., citet:king_2018_seasonal).

Ice velocity data are available with increasing spatial and temporal resolution (e.g., citet:vijay_2019_resolving). Until recently, ice velocity mosaics were limited to once per year during winter citep:joughin_2010_greenland, and they are still temporally limited, often to annual resolution, prior to 2000 (e.g. citet:mouginot_2018_1972to1990,mouginot_2018_1991to2000). Focusing on recent times, ice-sheet-wide velocity mosaics from the Sentinel-1A & 1B are now available every twelve days (http://PROMICE.org). The increased availability of satellite data has improved ice velocity maps both spatially and temporally, thereby decreasing the need to rely on spatial and temporal interpolation of velocities from annual/winter mosaics citep:andersen_2015_basin-scale,king_2018_seasonal,mouginot_2019_forty.

The discharge gates in this study are generated using only surface speed and an ice mask. We use the MEaSUREs Greenland Ice Sheet Velocity Map from InSAR Data, Version 2 citep:joughin_2010_greenland,NSIDC_0478, hereafter termed “MEaSUREs 0478” due to the National Snow and Ice Data Center (NSIDC) dateset ID number. We use the BedMachine v5 citep:morlighem_2017_bedmachine,NSIDC_BedMachine ice mask.

For ice thickness estimates, we use bed elevations from BedMachine v5 and the annual PROMICE PRODEM surface elevation (citet:winstrup_2024,winstrup_2024_data) from 2019 through 2022 (or the latest available). Prior to 2019 we adjust with surface elevation change from citet:khan_2016_geodetic. Ice sector and region delineation is from citet:mouginot_2019_glacier. Ice velocity data are obtained from a variety of products including Sentinel-1A & 1B derived by PROMICE citep:solgaard_2021, MEaSUREs 0478, MEaSUREs 0646 citep:NSIDC_0646, citet:mouginot_2018_1972to1990, and citet:mouginot_2018_1991to2000. Official glacier names come from citet:bjork_2015_brief. Other glacier names come from citet:mouginot_2019_glacier. See Table tab:data for an overview of datasets used in this work.

echo "times  all: " $(head -n1 ./tmp/dat_100_5000.csv | tr ',' '\n' | grep "vel_eff" | wc -l)
echo "times 19XX: " $(head -n1 ./tmp/dat_100_5000.csv | tr ',' '\n' | grep "vel_eff_19" | wc -l)
echo "times 20XX: " $(head -n1 ./tmp/dat_100_5000.csv | tr ',' '\n' | grep "vel_eff_20" |wc -l)
for Y in $(seq 2000 2025); do 
  echo "times ${Y}: " $(head -n1 ./tmp/dat_100_5000.csv | tr ',' '\n' | grep "vel_eff_${Y}" |wc -l)
done

This work uses src_bash[:eval yes]{head -n1 ./tmp/dat_100_5000.csv | tr ‘,’ ‘\n’ | grep “vel_eff” | wc -l} {{{results(3369)}}} different velocity maps, biased toward post-2015 when twelve-day ice velocities become available from the Sentinel-1 satellites. The temporal distribution is ~10 maps per year from 1986 to 2013, 14 in 2014, 25 in 2015, 38 in 2016, 81 in 2017, 54 in 2018, and one every ~12 days from 2019 onward.

Property	Name used in this paper	Reference
Basal topography	BedMachine v5	citet:morlighem_2017_bedmachine,NSIDC_BedMachine
Surface elevation	PROMICE PRODEM	citet:winstrup_2024,winstrup_2024_data
Surface elevation Change	PROMICE PRODEM 2019 – recent	citet:winstrup_2024,winstrup_2024_data
Surface elevation Change	Surface elevation change prior to 2019	citet:khan_2016_geodetic,GEUS_discharge_paper_elevation_change
Baseline velocity	MEaSUREs 0478	citet:NSIDC_0478
Velocity	PROMICE Sentinel	citet:solgaard_2021,solgaard_2021_data updated to 200 m
Velocity	MEaSUREs 0646	citet:NSIDC_0646
Velocity	MEaSUREs 0731	citet:NSIDC_0731,joughin_2010_greenland,joughin_2018_greenland
Velocity	pre-2000	citet:mouginot_2018_1972to1990,mouginot_2018_1991to2000
Sectors and regions	Sectors and regions	citet:mouginot_2019_glacier
Names		citet:bjork_2015_brief,mouginot_2019_glacier
Additional metadata		citet:Moon_2008,NSIDC_0642,Zwally_2012

Methods

Terminology

We use the following terminology, displayed in Fig. fig:overview:

“Pixels” are individual 200 m x 200 m raster discharge grid cells. We use the nearest neighbor when combining datasets that have different grid properties.
“Gates” are contiguous (including diagonal) clusters of pixels.
“Sectors” are spatial areas that have 0, 1, or > 1 gate(s) plus any upstream source of ice that flows through the gate(s), and come from citet:mouginot_2019_glacier.
“Regions” are groups of sectors, also from citet:mouginot_2019_glacier, and are labeled by approximate geographic region.
The “baseline” period is the average 2015, 2016, and 2017 winter velocity from MEaSUREs 0478.
“Coverage” is the percentage of total, region, sector, or gate discharge observed at any given time. By definition coverage is 100 % during the baseline period. From the baseline data, the contribution to total discharge of each pixel is calculated, and coverage is reported for all other maps that have missing observations (Fig. fig:coverage_schematic). Total estimated discharge is always reported because missing pixels are gap filled (see “Missing or invalid data” section below).
“Fast-flowing ice” is defined as ice that flows more than 100 m yr^-1.
Names are reported using the official Greenlandic names from citet:bjork_2015_brief; if an alternate name exists (e.g. from citet:mouginot_2019_glacier, or an English version), then this is shown in parentheses.

Although we refer to solid ice discharge, and it is in the solid phase when it passes the gates and eventually reaches the termini, submarine melting occurs at the termini and some of the discharge enters the fjord as liquid water citep:enderlin_2013_submarine.

Gate location

Gates are algorithmically generated for fast-flowing ice (greater than 100 m yr^-1) close to the ice sheet terminus determined by the baseline-period data. We apply a 2D inclusive mask to the baseline data for all ice flowing faster than 100 m yr^-1. We then select the mask edge where it is near the BedMachine ice mask (not including ice shelves), which effectively provides grounding line termini. We buffer the termini 5000 m in all directions creating ovals around the termini and once again down-select to fast-flowing ice pixels. This procedure results in gates 5000 m upstream from the baseline terminus that bisect the baseline fast-flowing ice. We manually mask some land- or lake-terminating glaciers which are initially selected by the algorithm due to fast flow and mask issues.

We select a 100 m yr^-1 speed cutoff because slower ice, taking longer to reach the terminus, is more influenced by SMB. For the influence of this threshold on our results see the Discussion section and Fig. fig:heatmap.

We select gates at 5000 m upstream from the baseline termini except at Sermeq Kujalleq (Jakobshavn Isbræ), which means that gates are likely > 5000 m from the termini further back in the historical record citep:murray_2015_extensive,wood_2018_ocean-induced. The choice of a 5000 m buffer follows from the fact that it is near terminus and thus avoids the need for (minor) SMB corrections downstream, yet is not too close to the terminus where discharge results are sensitive to the choice of distance-to-terminus value (Fig. fig:heatmap), which may be indicative of bed (ice thickness) errors. At Sermeq Kujalleq the termini has retreated ~5 km, so we move the baseline termini inland so that the final gate location is still a few km upstream of the present-day termini.

Thickness

We derive thickness from surface and bed elevation. We use the PROMICE PRODEM 2020 surface elevation as our baseline surface assigning it a fixed date of 2020-08-01, and the BedMachine bed elevations. We adjust the surface through time from 2020-08-01 by linearly interpolating elevation changes to other PRODEM years. We adjust the surface backward from 2019-08-01 using the citet:khan_2016 surface elevation change product, that covers the period from 1995 to 2019. We assume no surface change for dates prior to 1995 and from the most recent PRODEM to the most recent velocity timestamp. Finally, from the fixed bed and temporally varying surface, we calculate the time-dependent ice thickness at each gate pixel.

Missing or invalid data

The baseline data provide velocity at all gate locations by definition, but individual nonbaseline velocity maps often have missing or invalid data. Also, thickness provided by BedMachine is clearly incorrect in some places (e.g. fast-flowing ice that is 10 m thick, Fig. fig:h_v_histogram). We define invalid data and fill in missing data as described below.

Invalid velocity

We flag invalid (outlier) velocities by treating each pixel as an individual time series, applying a 30-point rolling window, flagging values more than 2 standard deviations outside the mean, and repeating this filter three times. We also drop the 1972 to 1985 years from citet:mouginot_2018_1972to1990 because there is low coverage and extremely high variability when using our algorithm.

This outlier detection method appears to correctly flag outliers (see citet:mankoff_2019_ice, for unfiltered time series graphs) but likely also flags some true short-term velocity increases. The effect of this filter is a ~1% reduction in discharge most years but more in years with high discharge – a reduction of 3.2 % in 2013, 4.3 % in 2003, and more in the 1980s when the data are noisy. Any analysis using these data and focusing on individual glaciers or short-term changes (or lack thereof) should reevaluate the upstream data sources.

Missing velocity

We generate an ice speed time series by assigning the PROMICE, MEaSUREs 0478, MEaSUREs 0646, and pre-2000 products to their respective reported time stamps (even though these are time-span products) or to the middle of their time span when they cover a long period such as the annual maps from citet:mouginot_2018_1972to1990,mouginot_2018_1991to2000. We ignore that any individual velocity map or pixel has a time span and not a time stamp. Velocities are sampled only where there are gate pixels. Missing pixel velocities are linearly interpolated in time, except for missing data at the beginning of the time series which are back- and forward filled with the temporally nearest value for that pixel (Fig. fig:coverage_schematic). We do not spatially interpolate missing velocities because the spatial changes around a missing data point are most likely larger than the temporal changes. We visually represent the discharge contribution of directly observed pixels, termed coverage (Fig. fig:coverage_schematic) as time series graphs and opacity of dots and error bars in the figures. The figures only display data where coverage is $≥$ 50 %, but the provided data files include coverage from 0 to 100 %. Therefore, the gap-filled discharge contribution at any given time is equal to 100 minus the coverage. Discharge is always reported as estimated total discharge even when coverage is less than 100 %.

Invalid thickness

The thickness data appear to be incorrect in some locations. For example, many locations have fast-flowing ice but report ice thickness as 10 m or less (Fig. fig:h_v_histogram, left panel). We accept all ice thickness greater than 20 m and construct from this a thickness vs. log₁₀-speed relationship. For all ice thickness less than or equal to 20 m thick we adjust thickness based on this relationship (Fig. fig:h_v_histogram, right panel). We selected the 20 m thickness cutoff after visually inspecting the velocity distribution (Fig. fig:h_v_histogram, left panel). This thickness adjustment adds 20 Gt yr^-1 to our baseline-period discharge estimate with no adjustment. In the Appendix and Table tab:thick_treatments we discuss the discharge contribution of these adjusted pixels, and a comparison among this and other thickness adjustments.

Discharge

We calculate discharge per pixel using density (917 kg m^-3), filtered and filled ice speed, projection-corrected pixel width, and adjusted ice thickness derived from time-varying surface elevation and a fixed bed elevation (Eq. eq:Q). We assume that any change in surface elevation corresponds to a change in ice thickness and thereby neglect basal uplift, erosion, and melt, which combined are orders of magnitude less than surface melting (e.g., citet:cowton_2012_rapid,khan_2007_elastic). We also assume depth-averaged ice velocity is equal to the surface velocity.

We calculate discharge using the gate orthogonal velocity at each pixel and at each timestamp – all velocity estimates are gate-orthogonal at all times, regardless of gate position, orientation, or changing glacier velocity direction over time.

Annual averages are calculated by linearly interpolating to daily and then estimating annual. The difference between this method and averaging only the observed samples is ~3 % median (5 % average, and a maximum of 10 % when examining the entire ice sheet and all years in our data). It is occasionally larger at individual glaciers when a year has few widely spaced samples of highly variable velocity.

Discharge uncertainty

\label{sec:D_uncertainty}

A longer discussion related to our and others treatments of errors and uncertainty is in the Appendix, but here we describe how we estimate the uncertainty related to the ice discharge following a simplistic approach. This yields an uncertainty of the total ice discharge of approximately 10 % throughout the time series.

At each pixel we estimate the maximum discharge, $D_\mathrm{max}$, from

\begin{equation} D_\mathrm{max} = ρ \, (V + σ_V) \, (H + σ_H) \, W, \end{equation}

and minimum discharge, $D_\mathrm{min}$, from

\begin{equation} D_\mathrm{min} = ρ \, (V - σ_V) \, (H - σ_H) \, W, \end{equation}

where $ρ$ is ice density, $V$ is baseline velocity, $σ_V$ is baseline velocity error, $H$ is ice thickness, $σ_H$ is ice thickness error, and $W$ is the width at each pixel. Included in the thickness term is surface elevation change through time ($\mathrm{d}H/\mathrm{d}t$). When datasets do not come with error estimates we treat the error as 0.

We use $ρ = 917$ kg m^-3 because the gates are near the terminus in the ablation zone, and ice thickness estimates should not include snow or firn, although regionally ice density may be < 917 kg m^-3 due to crevasses. We ignore the velocity error $σ_V$ because the proportional thickness error ($σ_H/H$) is an order of magnitude larger than the proportional velocity error ($σ_V/V$) yet both contribute linearly to the discharge. $W$ is location dependent due to the errors between our working map projection (EPSG 3413) and a more accurate spheroid model of the earth surface. We adjust linear gate width by up to ~4% in the north and ~-2.5% in the south of Greenland (area errors are up to 8%). On a pixel-by-pixel basis we used the provided thickness uncertainty except where we modified the thickness (H < 20 m); we prescribe an uncertainty of 0.5 times the adjusted thickness. Subsequently, the uncertainty on individual glacier, sector, region, or ice sheet scale is obtained by summing, but not reducing by the square of the sums, the uncertainty related to each pixel. We are conservative with our thickness error estimates, by assuming the uncertainty range is from $D_\mathrm{min}$ to $D_\mathrm{max}$ and not reducing by the sum of squares of sectors or regions.

Results

Gates

Our gate placement algorithm generates src_bash[:eval yes]{tail -n +2 tmp/dat_100_5000.csv | wc -l} {{{results(5851)}}} pixels making up src_bash[:eval yes]{cut -d”|” -f3 ./tmp/dat/gates_gateID@gates_100_5000.bsv | tail -n +2 | sort | uniq | wc -l} {{{results(0)}}} gates, assigned to src_bash[:eval yes]{cut -d”|” -f3 ./tmp/dat/sectors@Mouginot_2019.bsv | tail -n +2 | sort -n | uniq | wc -l} {{{results(173)}}} ice sheet sectors from citet:mouginot_2019_glacier. Previous similar studies have used 260 gates citep:mouginot_2019_forty, 230 gates citep:king_2018_seasonal, and 178 gates citep:enderlin_2014_improved.

The widest gate (~47 km) is Sermersuaq (Humboldt Gletsjer) and both Ikertivaq and Sermeq Kujalleq (Jakobshavn Isbræ) are ~34 km wide. A total of 14 glaciers have gate lengths longer than 10 km. The minimum gate width is 3 pixels (600 m) by definition in the algorithm.

The average unadjusted thickness gates is 399 m with a standard deviation of 253. The average thickness after adjustment is 429 m with a standard deviation of 223. A histogram of unadjusted and adjusted thickness at all gate locations is shown in Fig. fig:h_v_histogram.

Discharge

Our ice discharge dataset (Fig. fig:discharge_ts) reports a total discharge of 450 $±$ 43 Gt in 1986, has a minimum of 419 $±$ 39 Gt in 1996, and increases to 433 $±$ 40 in 2000 and further to 488 $±$ 45 Gt/yr in 2005, after which annual discharge remains approximately steady at 472 to 491 $±$ ~45 Gt/yr during the 2005 through 2020 period.

At the region scale, the SE glaciers (see Fig. fig:overview for regions) are responsible for 128 to 154 ($±$ 13 %) Gt yr^-1 of discharge (approximately one-third of ice-sheet-wide discharge) over the 1986 through 2019 period. By comparison, the predominantly land-terminating NO, NE, and SW together were also responsible for about one-third of total ice sheet discharge during this time (Fig. fig:discharge_ts_regions). The discharge from most regions has been approximately steady for the past decade. The NW region exhibited a persistent long-term increase in discharge – from ~90 to 114 Gt yr^-1 (27 % increase) over the 1999 through 2017 period, but has become more variable with declines and increases from 2017 through 2021. Increased variability also appears in the CW and CE regions beginning in 2016.

Focusing on eight major contributors at the individual sector or glacier scale (Fig. fig:discharge_ts_topfew), Sermeq Kujalleq (Jakobshavn Isbræ) has slowed down from an annual average high of ~50 Gt yr^-1 in 2013 to ~30 Gt yr^-1 in 2018, likely due to ocean cooling citep:khazendar_2019_interruption, but in 2021 returned briefly to nearly 50 Gt yr^-1. The Sermeq Kujalleq increasing trend and regular annual cycle has become disrupted in ~2015 with large decreases and shifting of the normal summer velocity maximum. Helheim briefly contributed more to sea level rise than Jakobshavn Isbræ in 2019, but has returned to 2nd place in 2020 and 2021 as Jakobshavn Isbræ speeds back up (Fig. fig:discharge_ts_topfew). We exclude Ikertivaq from the top eight because that gate spans multiple sectors and outlets, while the other top dischargers are each a single outlet.

Discussion

Different ice discharge estimates among studies likely stem from three categories: 1) changes in true discharge, 2) different input data (ice thickness and velocity), and 3) different assumptions and methods used to analyze data. Improved estimates of true discharge are the goal of this and many other studies, but changes in true discharge (category 1) can happen only when a work extends a time series into the future because historical discharge is fixed. Thus, any interstudy discrepancies in historical discharge must be due to category 2 (different data) or category 3 (different methods). Most studies use both updated data and new or different methods, but do not always provide sufficient information to disentangle the two. This is inefficient. To more quantitatively discuss interstudy discrepancies, it is imperative to explicitly consider all three potential causes of discrepancy. Only when results are fully reproducible – meaning all necessary data and code are available (cf. citet:mankoff_2017_past,rezvanbehbahani_2017_predicting,mankoff_2019_ice) – can new works confidently attribute discrepancies relative to old works. Therefore, in addition to providing new discharge estimates, we attempt to examine discrepancies among our estimates and other recent estimates. Without access to code and data from previous studies, it is challenging to take this examination beyond a qualitative discussion.

The algorithm-generated gates we present offer some advantages over traditional handpicked gates. Our gates are shared publicly, are generated by code that can be audited by others, and are easily adjustable within the algorithmic parameter space. This both allows sensitivity testing of gate location (Fig. fig:heatmap) and allows gate positions to systematically evolve with glacier termini (not done here).

Limitations

The discharge gates used here are generally ~5 km upstream of the glacier terminus. As glaciers retreat, they become closer to the terminus. For some glaciers that have retreated past the original gate location, the gates have been manually moved farther inland, meaning they are >5 km from the historical terminus. That ~5 km for the fastest flowing ice reaches the terminus within ~1/3 of a year. Meanwhile, the slowest flowing ice (100 m/yr cutoff, may slow down, usually speeds up) could take up to 50 years to reach the terminus. This means there is both a spatial and temporal lag between the observations reported here and the glacier terminus.

Due to this lag, this product is not suitable for processes studies at glacier termini that require high temporal resolution, unless it is combined with a terminus retreat product (e.g., citet:greene_2024). Recent work shows that this was negligible on an annual scale in the early part of this record (e.g., no net change, although there may still have been significant seasonal retreat and advance), but recently contributes an additional ~50 Gt yr^-1 (~10 %) mass loss from 2000 through 2020 citet:kochtitzky_2023,greene_2024.

We do not computer the SMB changes between the gate and the terminus due to the additional complexity required to do so, combined with their relatively small impact on the total discharge. Nonetheless, we note that this would contribute to double-counting some mass loss, if this discharge product is included unmodified in SMB mass loss estimates that include the entire ice sheet (also downstream of the gates used here). Recent work has quantified this effect, and we now report those values: It is estimated to be < 20 Gt yr^-1 by citet:kochtitzky_2022_mass (~4 % of ~500 Gt yr^-1 discharge), but we note citet:kochtitzky_2022_mass used flux gates closer to the terminus, meaning SMB losses in the product presented here would be larger.

Finally, it should be explicitly noted that what we term ice discharge here is divided into iceberg calving and submarine melt at the ice/ocean boundary. That division is highly variable in space and time, and without detailed work quantifying individual glaciers at specific points in time we estimate the split at ~50 % ± 40 % citep:enderlin_2013.

Data availability

This work in its entirety is available at doi:10.22008/promice/data/ice_discharge citep:GEUS_discharge_paper. The glacier-scale, sector, region, and Greenland summed ice sheet discharge dataset is available at doi:10.22008/promice/data/ice_discharge/d/v02 citep:GEUS_discharge_paper_d, where it will be updated as more velocity data become available. The gates can be found at doi:10.22008/promice/data/ice_discharge/gates/v02 citep:GEUS_discharge_paper_gates, the code at doi:10.22008/promice/data/ice_discharge/code/v0.0.1 citep:GEUS_discharge_paper_code, and the surface elevation change at doi:10.22008/promice/data/DTU/surface_elevation_change/v1.0.0 citep:GEUS_discharge_paper_elevation_change.

Conclusions

We have presented a novel dataset of flux gates and a 1986 through 2019 glacier-scale ice discharge estimate for the Greenland Ice Sheet. These data are underpinned by an algorithm that both selects gates for ice flux and then computes ice discharges.

Our results are similar to the most recent discharge estimate citep:king_2018_seasonal but begin in 1986 - although there are fewer samples prior to 2000. From our discharge estimate we show that over the past ~30 years, ice sheet discharge was ~440 Gt yr^-1 prior to 2000, rose to over 500 Gt yr^-1 from 2000 to 2005, and has held roughly steady since 2005 at near 500 Gt yr^-1. However, when viewed at a region or sector scale, the system appears more dynamic with spatial and temporal increases and decreases canceling each other out to produce the more stable ice sheet discharge. We note that there does not appear to be any dynamic connection among the regions, and any increase in one region that was offset by a decrease in another has likely been due to chance. If in coming years when changes occur the signals have matching signs, then ice sheet discharge would decrease or increase, rather than remain fairly steady.

The application of our flux gate algorithm shows that ice-sheet-wide discharge varies by ~30 Gt yr^-1 due only to gate position, or ~40 Gt yr^-1 due to gate position and cutoff velocity (Fig. fig:heatmap). This variance is approximately equal to the uncertainty associated with ice sheet wide discharge estimates reported in many studies (e.g. citet:rignot_2008_mass,andersen_2015_basin-scale,kjeldsen_2015_spatial). We highlight a major discrepancy with the ice discharge data of citet:enderlin_2014_improved and we suspect this discharge discrepancy – most pronounced in southeast Greenland – is associated with the choice of digital bed elevation model, specifically a deep hole in the bed at Køge Bugt.

Transparency in data and methodology are critical to move beyond a focus of estimating discharge quantities, towards more operational mass loss products with realistic errors and uncertainty estimates. The convention of devoting a paragraph, or even page, to methods is insufficient given the complexity, pace, and importance of Greenland Ice Sheet research citep:catania_2020. Therefore the flux gates, discharge data, and the algorithm used to generate the gates, discharge, and all figures from this paper are available. We hope that the flux gates, data, and code we provide here is a step toward helping others both improve their work and discover the errors in ours.

Other

References

Figures

Overview

Heatmap

Ice Thickness v. Velocity 2D Histogram: Color = count

Discharge Time Series

Discharge Time Series: Regions

Discharge Time Series: Top Few

Appendix

Errors and uncertainties

Here we describe our error and uncertainty treatments. We begin with a brief philosophical discussion of common uncertainty treatments, our general approach, and then the influence of various decisions made throughout our analysis, such as gate location and treatments of unknown thicknesses.

Traditional and mathematically valid uncertainty treatments divide errors into two classes: systematic (bias) and random. The primary distinction is that systematic errors do not decrease with more samples, and random errors decrease as the number of samples or measurements increases. The question is then which errors are systematic and which are random. A common treatment is to decide that errors within a region are systematic and among regions are random. This approach has no physical basis - two glaciers a few hundred meters apart but in different regions are assumed to have random errors, but two glaciers thousands of kilometers apart but within the same region are assumed to have systematic errors. It is more likely the case that all glaciers narrower than some width or deeper than some depth have systematic errors even if they are on opposite sides of the ice sheet, if ice thickness is estimated with the same method (i.e. the systematic error is likely caused by the sensor and airplane, not the location of the glacier).

The decision to have $R$ random samples (where $R$ is the number of regions, usually ~18 based on citet:zwally_2012_sectors) is also arbitrary. Mathematical treatment of random errors means that, even if the error is 50 % 18 measurements reduce it to only 11.79 %.

This reduction is unlikely to be physically meaningful. Our 173 sectors, 267 gates, and 5829 pixels means that, even if errors were 100 % for each, we could reduce it to 7.5, 6.1, or 1.3 % respectively. We note that the area error introduced by the common EPSG:3413 map projection is -5 % in the north and +8 % in the south. While this error is mentioned in some other works (e.g., citet:joughin_2018_greenland) it is often not explicitly mentioned.

We do not have a solution for the issues brought up here, except to discuss them explicitly and openly so that those, and our own, error treatments are clearly presented and understood to likely contain errors themselves.

Invalid thickness

We assume ice thicknesses < 20 m are incorrect where ice speed is > 100 m yr^-1. Of 5863 pixels, 5301 have valid thickness, and 562 (11 %) have invalid thickness. However, the speed at the locations of the invalid thicknesses is generally much less (and therefore the assumed thickness is less), and the influence on discharge is less than an average pixel with valid thickness (Table tab:thick_adjust).

When aggregating by gate, there are 267 gates. Of these, 181 (68 %) have no bad pixels and 86 (32 %) have some bad pixels, 61 have > 50 % bad pixels, and 57 (21 %) are all bad pixels.

We adjust these thickness using a poor fit (correlation coefficient: 0.3) of the log$₁₀$ of the ice speed to thickness where the relationship is known (thickness > 20 m). We set errors equal to one half the thickness (i.e. $σ_H = ± 0.5 \, H$). We also test the sensitivity of this treatment to simpler treatments, and have the following five categories:

NoAdj: No adjustments made. Assume BedMachine thicknesses are all correct.
300: If a gate has some valid pixel thicknesses, set the invalid thicknesses to the minimum of the valid thicknesses. If a gate has no valid thickness, set the thickness to 300 m.
400: Set all thicknesses < 50 m to 400 m
Fit: Use the thickness–speed relationship described above.

Table tab:thick_treatments shows the estimated baseline discharge to these four treatments:

Finally, Figure fig:gate_map shows the geospatial locations, concentration, and speed of gates with and without bad pixels.

Missing velocity

\label{sec:uncertain_vel}

We estimate discharge at all pixel locations for any time when there exists any velocity product. Not every velocity product provides velocity estimates at all locations, and we fill in where there are gaps by linearly interpolating velocity at each pixel in time. We calculate coverage, the discharge-weighted percent of observed velocity at any given time (Figure fig:coverage_schematic), and display coverage as 1) line plots over the time series graphs, 2) opacity of the error bars and 3) opacity of the infilling of time series dots. Linear interpolation and discharge-weighted coverage is illustrated in Figure fig:coverage_schematic, where pixel A has a velocity value at all three times, but pixel B has a filled gap at time $t_3$. The concentration of valid pixels is 0.5, but the weighted concentration, or coverage, is 9/11 or ~0.82. When displaying these three discharge values, $t_1$ and $t_4$ would have opacity of 1 (black), and $t_3$ would have opacity of 0.82 (dark gray).

This treatment is applied at the pixel level and then weight averaged to the gate, sector, region, and ice sheet results.

inkscape -z ./figs/gate_weight_schematic.svg -e ./figs/gate_weight_schematic.png

Errors from map projection

Our work takes place in a projected coordinate system (EPSG 3413) and therefore errors are introduced between the “true” earth spheroid (which is itself an approximation) and our projected coordinates system. We address these by calculating the projection error due to EPSG 3413 which is approximately +8 % in Northern Greenland and -6 % in Southern Greenland, and multiplying variables by a scaling factor if the variables do not already take this into account. Velocities are “true velocities” and not scaled, but the nominal 200 m gate width is scaled.

Velocity versus thickness

Software

This work was performed using only open-source software, primarily GRASS GIS citep:neteler_2012_GRASS and Python citep:van-rossum_1995_python, in particular the Jupyter citep:kluyver_2016_jupyter, pandas citep:mckinney_2010_pandas, numpy citep:oliphant_2006_numpy, statsmodel citep:seabold_2010_statsmodels, x-array citep:hoyer_2017_xarray, and Matplotlib citep:hunter_2007_matplotlib packages. The entire work was performed in Emacs citep:stallman_1981_emacs using Org Mode citep:schulte_2012_a-multi-language. The parallel citep:tange_2011_parallel tool was used to speed up processing. We used proj4 citep:proj4 to compute the errors in the EPSG 3413 projection. All code used in this work is available in the Supplemental Material.

Code

Misc Helper

Support pretty messages

RED='\033[0;31m'
ORANGE='\033[0;33m'
GREEN='\033[0;32m'
NC='\033[0m' # No Color
MSG_OK() { printf "${GREEN}${1}${NC}\n"; }
MSG_WARN() { printf "${ORANGE}WARNING: ${1}${NC}\n"; }
MSG_ERR() { echo "${RED}ERROR: ${1}${NC}\n" >&2; }

GRASS config

https://grass.osgeo.org/grass74/manuals/variables.html

GRASS_VERBOSE [all modules] toggles verbosity level -1 - complete silence (also errors and warnings are discarded) 0 - only errors and warnings are printed 1 - progress and important messages are printed (percent complete) 2 - all module messages are printed 3 - additional verbose messages are printed

export GRASS_VERBOSE=3
# export GRASS_MESSAGE_FORMAT=silent

if [ -z ${DATADIR+x} ]; then
    echo "DATADIR environment varible is unset."
    echo "Fix with: \"export DATADIR=/path/to/data\""
    exit 255
fi

set -x # print commands to STDOUT before running them

trap ctrl_c INT
function ctrl_c() {
  MSG_WARN "Caught CTRL-C"
  MSG_WARN "Killing process"
  kill -term $$ # send this program a terminate signal
}

Import Data

<<MSGS_pretty_print>>
<<GRASS_config>>

Bed and Surface

BedMachine v5

from Morlighem et al. (2017)

MSG_OK "BedMachine"
g.mapset -c BedMachine

for var in mask surface thickness bed errbed; do
  echo $var
  r.external source=NetCDF:${DATADIR}/Morlighem_2017/BedMachineGreenland-v5.nc:${var} output=${var}
done

r.colors -a map=errbed color=haxby

g.mapset PERMANENT
g.region raster=surface@BedMachine res=200 -a -p
g.region -s
g.mapset BedMachine
g.region -dp

r.colors map=mask color=haxby

r.mapcalc "mask_ice = if(mask == 2, 1, null())"

Sectors

Mouginot 2019

From citet:mouginot_2019_glacier

Import & Clean

MSG_OK "Mouginot 2019 sectors"

g.mapset -c Mouginot_2019
v.in.ogr input=${DATADIR}/Mouginot_2019 output=sectors_all
v.extract input=sectors_all where="NAME NOT LIKE '%ICE_CAP%'" output=sectors

db.select table=sectors | head
v.db.addcolumn map=sectors columns="region_name varchar(100)"
db.execute sql="UPDATE sectors SET region_name=SUBREGION1 || \"___\" || NAME"

v.to.db map=sectors option=area columns=area units=meters

mkdir -p ./tmp/
# db.select table=sectors > ./tmp/Mouginot_2019.txt

v.to.rast input=sectors output=sectors use=cat label_column=region_name
r.mapcalc "mask_GIC = if(sectors)"

# # regions map
v.to.rast input=sectors output=regions_tmp use=cat label_column=SUBREGION1
# which categories exist?
# r.category regions separator=comma | cut -d, -f2 | sort | uniq
# Convert categories to numbers
r.category regions_tmp separator=comma | sed s/NO/1/ | sed s/NE/2/ | sed s/CE/3/ | sed s/SE/4/ | sed s/SW/5/ | sed s/CW/6/ | sed s/NW/7/ > ./tmp/mouginot.cat
r.category regions_tmp separator=comma rules=./tmp/mouginot.cat
# r.category regions_tmp
r.mapcalc "regions = @regions_tmp"

# # region vector 
# r.to.vect input=regions output=regions type=area
# v.db.addcolumn map=regions column="REGION varchar(2)"
# v.what.vect map=regions column=REGION query_map=sectors query_column=SUBREGION1

# # mask

Test

grass74 ./G/Mouginot_2019
d.mon start=wx0
d.rast regions
d.rast sectors
d.vect sectors_all fill_color=none color=red
d.vect sectors fill_color=none

Zwally 2012

I use an “expanded boundary” version. This was created by loading the Zwally sectors into QGIS and moving the coasts outward. This is done because some gates (ice) is outside the boundaries provided by Zwally.

g.mapset -c Zwally_2012
v.in.ogr input=${DATADIR}/Zwally_2012/sectors_enlarged output=Zwally_2012

2D Area Error

EPSG:3413 has projection errors of $±$ ~8% in Greenland
Method
- Email: Re: {GRASS-user} scale error for each pixel
- Webmail: https://www.mail-archive.com/grass-user@lists.osgeo.org/msg35005.html

MSG_OK "2D Area Error"
g.mapset PERMANENT

if [[ "" == $(g.list type=raster pattern=err_2D) ]]; then
    r.mask -r
    g.region -d

    g.region res=1000 -ap # do things faster
    r.mapcalc "x = x()"
    r.mapcalc "y = y()"
    r.latlong input=x output=lat_low
    r.latlong -l input=x output=lon_low

    r.out.xyz input=lon_low,lat_low separator=space > ./tmp/llxy.txt
    PROJSTR=$(g.proj -j)
    echo $PROJSTR

    paste -d" " <(cut -d" " -f1,2 ./tmp/llxy.txt) <(cut -d" " -f3,4 ./tmp/llxy.txt | proj -VS ${PROJSTR} | grep Areal | column -t | sed s/\ \ /,/g | cut -d, -f4) > ./tmp/xy_err.txt

    r.in.xyz input=./tmp/xy_err.txt  output=err_2D_inv separator=space
    r.mapcalc "err_2D = 1/(err_2D_inv^0.5)" # convert area error to linear multiplier error
    g.region -d

    r.latlong input=x output=lat # for exporting at full res
    r.latlong -l input=x output=lon
fi

# sayav done
g.region -d

Velocity

MEaSUREs

See: MEaSUREs README

[X] 0478: 2000 – 2017 annual average
[X] 0481: 6-11 day velocity TSX
[X] 0646: Monthly velocity - sparse glacier coverage 1985 through 2016
[ ] 0670: 1995 – 2015 average
[ ] 0725: 2015 & 2016 annual average
[X] 0731: Monthly ice sheet velocity 2015 through 2018
[X] 0766: S1 6-d16

0478.002

MEaSUREs Greenland Ice Sheet Velocity Map from InSAR Data, Version 2
Winter velocity maps

Import

First read in the 200 m files
Then read in the 500 m files if there were no 200 m files

MSG_OK "MEaSURES.0478"
g.mapset -c MEaSUREs.0478

MSG_OK "  200 m..."
r.mask -r
ROOT=${DATADIR}/MEaSUREs/NSIDC-0478.002/
VX=$(find ${ROOT} -name "*mosaic200_*vx*.tif" | head -n1) # DEBUG
for VX in $(find ${ROOT} -name "*mosaic200_*vx*.tif" | LC_ALL=C sort); do
  VY=${VX/vx/vy}
  EX=${VX/vx/ex}
  EY=${EX/ex/ey}
  DATE=$(dirname ${VX} | rev | cut -d"/" -f1 | rev | sed s/\\./_/g)
  # echo $DATE
  # need to import not link to external so that we can set nulls to 0
  parallel --verbose --bar r.in.gdal input={1} output={2}_${DATE} ::: ${VX} ${VY} ${EX} ${EY} :::+ VX VY EX EY
  parallel --verbose --bar r.null map={}_${DATE} null=0 ::: VX VY EX EY
done
g.region raster=VX_${DATE} -pa

MSG_OK "  500 m..."
VX=$(find ${ROOT} -name "*mosaic500_*vx*.tif" | head -n1) # DEBUG
for VX in $(find ${ROOT} -name "*mosaic500_*vx*.tif" | LC_ALL=C sort); do
  VY=${VX/vx/vy}
  EX=${VX/vx/ex}
  EY=${EX/ex/ey}
  DATE=$(dirname ${VX} | rev | cut -d"/" -f1 | rev | sed s/\\./_/g)
  echo $DATE

  # Read in all the 500 m velocity data
  parallel --verbose --bar r.external source={1} output={2}_${DATE}_500 ::: ${VX} ${VY} ${EX} ${EY} :::+ VX VY EX EY 
  # If the 200 m data exists, will produce an error and continue
  # If the 200 m data does not exist, will resample from 500
  r.mapcalc "VX_${DATE} = VX_${DATE}_500"
  r.mapcalc "VY_${DATE} = VY_${DATE}_500"
  r.mapcalc "EX_${DATE} = EX_${DATE}_500"
  r.mapcalc "EY_${DATE} = EY_${DATE}_500"
  parallel --verbose --bar r.null map={}_${DATE} null=0 ::: VX VY EX EY
done

Baseline: Average of 2015-2017

See ./dat/remove_ice_manual.kml
This is due to extensive Jakobshavn retreat between baseline and present
The gates need to be >5 km from the baseline terminus

MSG_OK "Baseline"
g.mapset -c MEaSUREs.0478

r.series input=VX_2015_09_01,VX_2016_09_01,VX_2017_09_01 output=vx_baseline method=average range=-1000000,1000000
r.series input=VY_2015_09_01,VY_2016_09_01,VY_2017_09_01 output=vy_baseline method=average range=-1000000,1000000

r.series input=EX_2015_09_01,EX_2016_09_01,EX_2017_09_01 output=ex_baseline method=average range=-1000000,1000000
r.series input=EY_2015_09_01,EY_2016_09_01,EY_2017_09_01 output=ey_baseline method=average range=-1000000,1000000

v.import input=./dat/remove_ice_manual.kml output=remove_ice_manual --o
r.mask -i vector=remove_ice_manual --o

r.mapcalc "vel_baseline = sqrt(vx_baseline^2 + vy_baseline^2)"
r.mapcalc "vel_err_baseline = sqrt(ex_baseline^2 + ey_baseline^2)"

r.mask -r

parallel --verbose --bar r.null map={}_baseline setnull=0 ::: vx vy vel ex ey vel_err
r.colors -e map=vel_baseline,vel_err_baseline color=viridis

Fill in holes

There are holes in the velocity data which will create false gates. Fill them in.
Clump based on yes/no velocity
- Largest clump is GIS
- 2nd largest is ocean
Mask by ocean (so velocity w/ holes remains)
Fill holes

r.mask -r
r.mapcalc "no_vel = if(isnull(vel_baseline), 1, null())"
r.mask no_vel
r.clump input=no_vel output=no_vel_clump --o
ocean_clump=$(r.stats -c -n no_vel_clump sort=desc | head -n1 | cut -d" " -f1)
r.mask -i raster=no_vel_clump maskcats=${ocean_clump} --o
r.fillnulls input=vel_baseline out=vel_baseline_filled method=bilinear
r.mask -r
g.rename raster=vel_baseline_filled,vel_baseline --o
r.colors map=vel_baseline -e color=viridis

Display

d.mon start=wx0
d.erase
d.rast vel
d.rast vel_filled

0481.004 TSX

Generate VRTs

One map per date
Build GDAL virtual tiles for every date (when data exists)

g.mapset -c MEaSUREs.0481

ROOT=${DATADIR}/MEaSUREs/NSIDC-0481.004
VRTROOT=./tmp/NSIDC-0481.004.vrt/
mkdir -p ${VRTROOT}

for date in $(cd ${ROOT}; ls -d */|tr -d '/'); do
  if [[ ! -f ${VRTROOT}/${date}_vx.vrt ]]; then # VRT file does not exist?
    LIST=$(find ${ROOT}/${date}/ -name "TSX*vx*.tif" | LC_ALL=C sort)
    if [[ ! -z ${LIST} ]]; then
      MSG_OK "Building VRTs for ${date}"

      VX=$(ls $ROOT/$date/*vx* | head -n1)
      T0=$(basename ${VX} | cut -d_ -f3)
      T1=$(basename ${VX} | cut -d_ -f4)
      SEC0=$(date --utc --date="${T0}" +"%s")
      SEC1=$(date --utc --date="${T1}" +"%s")
      MID=$(echo "(${SEC0}+${SEC1})/2"|bc)
      DATE=$(date --utc --date="@${MID}" +"%Y_%m_%d")
      
      parallel --verbose --bar gdalbuildvrt -overwrite ${VRTROOT}/${DATE}_{}.vrt $\(find ${ROOT}/${date} -name "TSX*{}*.tif" \| LC_ALL=C sort\) ::: vx vy ex ey
    fi
  fi
done

Import VRTs

MSG_OK "MEaSURES.0481"
g.mapset -c MEaSUREs.0481

# Set a super region
gdalbuildvrt -overwrite tmp/0481.vrt ./tmp/NSIDC-0481.004.vrt/*_vx.vrt
r.external source=tmp/0481.vrt output=tmp -e -r 
g.region raster=tmp -pa res=250
g.remove -f type=raster name=tmp

r.mask -r
ROOT=./tmp/NSIDC-0481.004.vrt/
VX=$(find ${ROOT} -name "*vx*.vrt" | head -n1) # debug
for VX in $(find ${ROOT} -name "*vx*.vrt" | LC_ALL=C sort); do
    VY=${VX/vx/vy}
    EX=${VX/vx/ex}
    EY=${EX/ex/ey}
    DATE=$(basename $VX | cut -d"_" -f1-3)
    echo $DATE
    
    parallel --verbose --bar r.external source={1} output={2}_${DATE} ::: ${VX} ${VY} ${EX} ${EY} :::+ VX VY EX EY
done

0646.003

MEaSUREs Greenland Ice Velocity: Selected Glacier Site Velocity Maps from Optical Images, Version 2
Monthly velocity maps

Generate VRTs

One map per month
Build GDAL virtual tiles for every month (when data exists)

g.mapset -c MEaSUREs.0646

ROOT=${DATADIR}/MEaSUREs/NSIDC-0646.003/
VRTROOT=./tmp/NSIDC-0646.003.vrt/
mkdir -p ${VRTROOT}
for year in $(seq 1985 2018); do
  for month in $(seq -w 1 12); do
    if [[ ! -f ${VRTROOT}/${year}_${month}_vx.vrt ]]; then # VRT file does not exist?
      LIST=$(find ${ROOT} -name "*${year}-${month}_vx_*.tif" | LC_ALL=C sort)
      if [[ ! -z ${LIST} ]]; then
        MSG_OK "Building VRTs for ${year} ${month}"
        parallel --verbose --bar gdalbuildvrt -overwrite ${VRTROOT}/${year}_${month}_{}.vrt $\(find ${ROOT} -name "*${year}-${month}_{}_*.tif" \| LC_ALL=C sort\) ::: vx vy ex ey
      fi
    fi
  done
done

Import VRTs

MSG_OK "MEaSURES.0646"
g.mapset -c MEaSUREs.0646

r.mask -r
ROOT=./tmp/NSIDC-0646.003.vrt/
VX=$(find ${ROOT} -name "*vx*.vrt" | head -n1) # debug
for VX in $(find ${ROOT} -name "*vx*.vrt" | LC_ALL=C sort); do
    VY=${VX/vx/vy}
    EX=${VX/vx/ex}
    EY=${EX/ex/ey}
    DATE=$(basename $VX | cut -d"_" -f1-2)
    DATE=${DATE}_15
    echo $DATE
    
    parallel --verbose --bar r.external source={1} output={2}_${DATE} ::: ${VX} ${VY} ${EX} ${EY} :::+ VX VY EX EY
done
g.region raster=VX_${DATE} -pa
# g.list type=raster mapset=MEaSUREs.0646

0731.005

MEaSUREs Greenland Monthly Ice Sheet Velocity Mosaics from SAR and Landsat, Version 1

Import

MSG_OK "MEaSURES.0731"
g.mapset -c MEaSUREs.0731
r.mask -r
ROOT=${DATADIR}/MEaSUREs/NSIDC-0731.005/
VX=$(find ${ROOT} -name "*mosaic_*vx*.tif" | head -n1) # DEBUG
for VX in $(find ${ROOT} -name "*mosaic_*vx*.tif" | LC_ALL=C sort); do
  VY=${VX/vx/vy}
  EX=${VX/vx/ex}
  EY=${EX/ex/ey}

  T0=$(dirname ${VX} | rev | cut -d"/" -f1 | rev|cut -d"_" -f4 | tr '.' '-')
  T1=$(dirname ${VX} | rev | cut -d"/" -f1 | rev|cut -d"_" -f5 | tr '.' '-')
  SEC0=$(date --utc --date="${T0}" +"%s")
  SEC1=$(date --utc --date="${T1}" +"%s")
  MID=$(echo "(${SEC0}+${SEC1})/2"|bc)
  DATE=$(date --utc --date="@${MID}" +"%Y_%m_%d")

  # echo $DATE
  parallel --verbose --bar r.external source={1} output={2}_${DATE} ::: ${VX} ${VY} ${EX} ${EY} :::+ VX VY EX EY
  parallel --verbose --bar r.null map={}_${DATE} null=0 ::: VX VY EX EY
done
g.region raster=VX_${DATE} -pa

0766.002 Sentinel

Import

MSG_OK "MEaSURES.0766"
g.mapset -c MEaSUREs.0766
r.mask -r
ROOT=${DATADIR}/MEaSUREs/NSIDC-0766.002
VX=$(find ${ROOT} ${ROOT}-updates -name "*mosaic_*vx*.tif" | tail -n1) # DEBUG
for VX in $(find ${ROOT} ${ROOT}-updates -name "*mosaic_*vx*.tif" | LC_ALL=C sort); do
  VY=${VX/vx/vy}
  EX=${VX/vx/ex}
  EY=${EX/ex/ey}

  T0=$(basename ${VX} | cut -d_ -f5)
  T1=$(basename ${VX} | cut -d_ -f6)
  SEC0=$(date --utc --date="${T0}" +"%s")
  SEC1=$(date --utc --date="${T1}" +"%s")
  MID=$(echo "(${SEC0}+${SEC1})/2"|bc)
  DATE=$(date --utc --date="@${MID}" +"%Y_%m_%d")

  # echo $DATE
  parallel --verbose --bar r.external source={1} output={2}_${DATE} ::: ${VX} ${VY} ${EX} ${EY} :::+ VX VY EX EY
  # Can't r.null for external maps.
  # parallel --verbose --bar r.null map={}_${DATE} null=0 ::: VX VY EX EY 
done
g.region raster=VX_${DATE} -pa

PROMICE IV 200m

Data Intro

DIR=${DATADIR}/Promice200m/
(cd ${DIR}; ls *.nc | head)
(cd ${DIR}; ncdump -h $(ls *.nc | head -n1) | grep "float")

Import data

Read in all the data
Conversion from [m day-1] to [m year-1] is done in section Just one velocity cutoff & buffer distance

MSG_OK "Promice 200m"
g.mapset -c promice
ROOT=${DATADIR}/Promice200m/

FILE=$(find ${ROOT} -name "*.nc" | head -n1) # DEBUG 
for FILE in $(find ${ROOT} -name "*.nc" | LC_ALL=C sort); do
  DATE_STR=$(ncdump -t -v time ${FILE} | tail -n2 | tr -dc '[0-9\-]' | tr '-' '_')
  echo $DATE_STR

  r.external -o source="NetCDF:${FILE}:land_ice_surface_easting_velocity" output=vx_${DATE_STR}
  r.external -o source="NetCDF:${FILE}:land_ice_surface_northing_velocity" output=vy_${DATE_STR}

  r.external -o source="NetCDF:${FILE}:land_ice_surface_easting_velocity_std" output=ex_${DATE_STR}
  r.external -o source="NetCDF:${FILE}:land_ice_surface_northing_velocity_std" output=ey_${DATE_STR}
done

Mouginot 2018 (pre-2000 velocities)

See citet:mouginot_2018_1972to1990 and citet:mouginot_2018_1991to2000

MSG_OK "Mouginot pre 2000"
g.mapset -c Mouginot_pre2000

ROOT=${DATADIR}/Mouginot_2018/D1GW91
find ${ROOT} -name "*.nc"
FILE=$(find ${ROOT} -name "*.nc" | head -n1 | LC_ALL=C sort) # DEBUG
for FILE in $(find ${ROOT} -name "*.nc"); do
  YYYYMMDD=$(echo ${FILE} | cut -d"_" -f4)
  YEAR=$(echo ${YYYYMMDD} | cut -d"-" -f1)
  DATE=${YEAR}_01_01
  echo $DATE
  r.external -o source="NetCDF:${FILE}:VX" output=vx_${DATE}
  r.external -o source="NetCDF:${FILE}:VY" output=vy_${DATE}
done

# ROOT=${DATADIR}/Mouginot_2018/D1MM37
# find ${ROOT} -name "*.nc"
# FILE=$(find ${ROOT} -name "*.nc" | head -n1) # DEBUG
# for FILE in $(find ${ROOT} -name "*.nc"); do
#   YYYYMMDD=$(echo ${FILE} | cut -d"_" -f4)
#   YEAR=$(echo ${YYYYMMDD} | cut -d"-" -f1)
#   DATE=${YEAR}_01_01
#   echo $DATE
#   r.external -o source="NetCDF:${FILE}:VX" output=vx_${DATE}
#   r.external -o source="NetCDF:${FILE}:VY" output=vy_${DATE}
# done

Display

d.mon start=wx0
g.list type=raster pattern=vx_*

d.erase; d.rast vx_1990-07-01
d.erase; d.rast vx_1991-07-01
d.erase; d.rast vx_1992-07-01
d.erase; d.rast vx_1993-07-01
d.erase; d.rast vx_1994-07-01
d.erase; d.rast vx_1995-07-01
d.erase; d.rast vx_1996-07-01
d.erase; d.rast vx_1997-07-01
d.erase; d.rast vx_1998-07-01
d.erase; d.rast vx_1999-07-01

Glacier Names

From Bjørk et al. (2015).
Also use citet:mouginot_2019_glacier

Bjørk 2015

Write out x,y,name. Can use x,y and mean gate location to find closest name for each gate.

MSG_OK "Bjørk 2015"
g.mapset -c Bjork_2015

ROOT=${DATADIR}/Bjørk_2015/

cat ${ROOT}/GreenlandGlacierNames_GGNv01.csv |  iconv -c -f utf-8 -t ascii | grep GrIS | awk -F';' '{print $3"|"$2"|"$7}' | sed s/,/./g | m.proj -i input=- | sed s/0.00\ //g | v.in.ascii input=- output=names columns="x double precision, y double precision, name varchar(99)"

# db.select table=names | tr '|' ',' > ./tmp/Bjork_2015_names.csv

Mouginot 2019

g.mapset Mouginot_2019
db.select table=sectors | head
# v.out.ascii -c input=sectors output=./tmp/Mouginot_2019_names.csv columns=NAME,SUBREGION1

NSIDC 0642

https://nsidc.org/data/NSIDC-0642/versions/1

MSG_OK "NSIDC 0642"
g.mapset -c NSIDC_0642
ROOT=${DATADIR}/Moon_2008
v.in.ogr input=${ROOT}/GlacierIDs_v01.2.shp output=GlacierIDs
# v.db.select map=GlacierIDs | head

ROOT=${DATADIR}/Moon_2008
md5sum ${ROOT}/GlacierIDs_v01.2.*

Elevation

h_0 is PRODEM 20 set to day 182 of 2020
h_n+ is PRODEM following years
- When PRODEM ends, continue using last year (e.g., pandas ffill())
h_n- is Khan 2016e ${DATADIR}/Khan_2016/README.org
- Coverage is 1995 through 2020
- Prior to Khan beginning, use first year (e.g., pandas bfill())

PRODEM

MSG_OK "dh/dt"

g.mapset -c PRODEM
r.mask -r

f=$(ls ${DATADIR}/PRODEM/PRODEM??.tif | head -n1) # debug
for f in $(ls ${DATADIR}/PRODEM/PRODEM??.tif); do
  y=20$(echo ${f: -6:2})
  r.in.gdal -r input=${f} output=DEM_${y} band=1
  # r.in.gdal -r input=${f} output=var_${y} band=2
  # r.in.gdal -r input=${f} output=dh_${y} band=3
  # r.in.gdal -r input=${f} output=time_${y} band=4
  # r.univar -g time_2019 # mean = DOI 213 = 01 Aug
done
g.region raster=DEM_2019 -pa

Khan 2016

MSG_OK "Khan 2016"

g.mapset -c Khan_2016
r.mask -r

g.region -d
g.region res=2000 -pa

cat << EOF > ./tmp/elev_filter.txt
TITLE     See r.mfilter manual
    MATRIX    3
    1 1 1
    1 1 1
    1 1 1
    DIVISOR   0
    TYPE      P
EOF

FILE=${DATADIR}/Khan_2016/dhdt_1995-2015_GrIS.txt
head -n1 $FILE
Y=1995 # debug
for Y in $(seq 1995 2010); do
  col=$(echo "$Y-1995+3" | bc -l)
  echo $Y $col
  if [[ "" == $(g.list type=raster pattern=dh_${Y}) ]]; then
    # remove comments, leading spaces, and convert
    # spaces to comma, swap lat,lon, then import
    cat ${FILE} \
      | grep -v "^%" \
      | sed s/^\ *//g \
      | sed s/\ \ \*/,/g \
      | cut -d"," -f1,2,${col} \
      | awk -F, '{print $2 "|" $1 "|" $3}' \
      | m.proj -i input=- \
      | r.in.xyz input=- output=dh_${Y}_unfiltered
  fi
done


FILE=${DATADIR}/Khan_2016/GR_2011_2020.txt
head -n6 $FILE
Y=2011 # debug
for Y in $(seq 2011 2019); do
  col=$(echo "($Y-2011)*2 +3" | bc -l)
  echo $Y $col
  if [[ "" == $(g.list type=raster pattern=dh_${Y}) ]]; then
    # remove comments, leading spaces, and convert
    # spaces to comma, swap lat,lon, then import
    cat ${FILE} \
      | grep -v "^%" \
      | sed s/^\ *//g \
      | sed s/\ \ \*/,/g \
      | cut -d"," -f1,2,${col} \
      | awk -F, '{print $2 "|" $1 "|" $3}' \
      | m.proj -i input=- \
      | r.in.xyz input=- output=dh_${Y}_unfiltered
  fi
done

parallel "r.mfilter -z input=dh_{1}_unfiltered output=dh_{1} filter=./tmp/elev_filter.txt repeat=2" ::: $(seq 1995 2019)
parallel "r.colors map=dh_{1} color=difference" ::: $(seq 1995 2019)

DEM

Merge Khan dh/dt w/ PRODEM to generate annual DEMs

MSG_OK "DEM"
g.mapset -c DEM

g.region raster=DEM_2020@PRODEM -pa
for y in {2019..2023}; do
  r.mapcalc "DEM_${y} = DEM_${y}@PRODEM"
done

for y in {2019..1995}; do
  y1=$(( ${y} + 1 ))
  r.mapcalc "DEM_${y} = DEM_${y1} - dh_${y}@Khan_2016"
done

Find Gates

Algorithm

[X] Find all fast-moving ice (>X m yr^-1)
- Results not very sensitive to velocity limit (10 to 100 m yr^-1 examined)
[X] Find grounding line by finding edge cells where fast-moving ice borders water or ice shelf based (loosely) on BedMachine mask
[X] Move grounding line cells inland by X km, again limiting to regions of fast ice.
- Results not very sensitive to gate position (1 - 5 km range examined)
[X] Discard gates if group size ∈ [1,2]
[X] Manually clean a few areas (e.g. land-terminating glaciers, gates due to invalid masks, etc.) by manually selecting invalid regions in Google Earth, then remove gates in these regions

Note that “fast ice” refers to flow velocity, not the sea ice term of “stuck to the land”.

INSTRUCTIONS: Set VELOCITY_CUTOFF and BUFFER_DIST to 50 and 2500 respectively and run the code. Then repeat for a range of other velocity cutoffs and buffer distances to get a range of sensitivities.

OR: Tangle via ((org-babel-tangle) the code below (C-c C-v C-t or ) to ./gate_IO.sh and then run this in a GRASS session:

<<MSGS_pretty_print>>
<<GRASS_config>>

VELOCITY_CUTOFF=100
BUFFER_DIST=5000
. ./gate_IO.sh

Create a new mapset for this specific velocity cutoff and buffer distance

g.mapset -c gates_vel_buf
g.region -d

From above:

[X] Find grounding line by finding edge cells where fast-moving ice borders water or ice shelf based (loosely) on BedMachine mask

The “loosely” is because the BedMachine mask doesn’t always reach into each fjord all the way. I buffer the BedMachine mask by 2 km here so that it extends to the edge of the velocity data.

g.copy raster=mask_ice@BedMachine,mask_ice --o
# Grow by 2 km (10 cells @ 200 m/cell)
r.grow input=mask_ice output=mask_ice_grow radius=10 new=1 --o
r.mask mask_ice_grow

The fast ice edge is where there is fast-flowing ice overlapping with not-ice.

r.mapcalc "fast_ice = if(vel_baseline@MEaSUREs.0478 > ${VELOCITY_CUTOFF}, 1, null())" --o
r.mask -r

# no velocity data, or is flagged as ice shelf or land in BedMachine
r.mapcalc "not_ice = if(isnull(vel_baseline@MEaSUREs.0478) ||| (mask@BedMachine == 0) ||| (mask@BedMachine == 3), 1, null())" --o

r.grow input=not_ice output=not_ice_grow radius=1.5 new=99 --o
r.mapcalc "fast_ice_edge = if(((not_ice_grow == 99) && (fast_ice == 1)), 1, null())" --o

The gates are set ${BUFFER_DIST} inland from the fast ice edge. This is done by buffering the fast ice edge (which fills the space between the fast ice edge and buffer extent) and then growing the buffer by 1. This last step defines the gate locations.

However, in order to properly estimate discharge, the gate location is not enough. Ice must flow from outside the gates, through the gates, to inside the gates, and not flow from one gate pixel to another gate pixel (or it would be counted 2x).

r.buffer input=fast_ice_edge output=fast_ice_buffer distances=${BUFFER_DIST} --o
r.grow input=fast_ice_buffer output=fast_ice_buffer_grow radius=1.5 new=99 --o
r.mask -i not_ice --o
r.mapcalc "gates_inside = if(((fast_ice_buffer_grow == 99) && (fast_ice == 1)), 1, null())" --o
r.mask -r

r.grow input=gates_inside output=gates_inside_grow radius=1.1 new=99 --o
r.mask -i not_ice --o
r.mapcalc "gates_maybe = if(((gates_inside_grow == 99) && (fast_ice == 1) && isnull(fast_ice_buffer)), 1, null())" --o
r.mask -r

r.grow input=gates_maybe output=gates_maybe_grow radius=1.1 new=99 --o
r.mask -i not_ice --o
r.mapcalc "gates_outside = if(((gates_maybe_grow == 99) && (fast_ice == 1) && isnull(fast_ice_buffer) && isnull(gates_inside)), 1, null())" --o
r.mask -r

r.mapcalc "gates_IO = 0" --o
r.mapcalc "gates_IO = if(isnull(gates_inside), gates_IO, 1)" --o
r.mapcalc "gates_IO = if(isnull(gates_outside), gates_IO, -1)" --o

r.colors map=gates_inside color=red
r.colors map=gates_maybe color=grey
r.colors map=gates_outside color=blue
r.colors map=gates_IO color=viridis

For each gate, split into two for the vector components of the velocity, then…
If flow is from gate to INSIDE, it is discharged
If flow is from gate to GATE, it is ignored
If flow is from gate to NOT(GATE || INSIDE) it is ignored
- If gates are a closed loop, such as the 1700 m flight-line, then this scenario would be NEGATIVE discharge, not ignored. This was tested with the 1700 m flight-line and compared against both the vector calculations and WIC estimates.

var	value	meaning
vx	> 0	east / right
vx	< 0	west / left
vy	> 0	north / up
vy	< 0	south / down
GRASS indexing	[0,1]	cell to the right
	[0,-1]	left
	[-1,0]	above
	[1,0]	below

# g.mapset -c gates_50_2500

r.mask -r

r.mapcalc "gates_x = 0" --o
r.mapcalc "gates_x = if((gates_maybe == 1) && (vx_baseline@MEaSUREs.0478 > 0), gates_IO[0,1], gates_x)" --o
r.mapcalc "gates_x = if((gates_maybe != 0) && (vx_baseline@MEaSUREs.0478 < 0), gates_IO[0,-1], gates_x)" --o

r.mapcalc "gates_y = 0" --o
r.mapcalc "gates_y = if((gates_maybe != 0) && (vy_baseline@MEaSUREs.0478 > 0), gates_IO[-1,0], gates_y)" --o
r.mapcalc "gates_y = if((gates_maybe != 0) && (vy_baseline@MEaSUREs.0478 < 0), gates_IO[1,0], gates_y)" --o

r.mapcalc "gates_x = if(gates_x == 1, 1, 0)" --o
r.mapcalc "gates_y = if(gates_y == 1, 1, 0)" --o

r.null map=gates_x null=0 # OR r.null map=gates_x setnull=0
r.null map=gates_y null=0 # OR r.null map=gates_y setnull=0

Clean Gates

Subset to where there is known discharge

r.mapcalc "gates_xy_clean00 = if((gates_x == 1) || (gates_y == 1), 1, null())" --o
r.mapcalc "gates_xy_clean0 = if(gates_xy_clean00 & if(DEM_2019@DEM), 1, null())" --o

Remove small areas (clusters <X cells)

# Remove clusters of 2 or less. How many hectares in X pixels?
# frink "(200 m)^2 * 2 -> hectares" # ans: 8.0

r.clump -d input=gates_xy_clean0 output=gates_gateID --o
r.reclass.area -d input=gates_gateID output=gates_area value=9 mode=lesser method=reclass --o

r.mapcalc "gates_xy_clean1 = if(isnull(gates_area), gates_xy_clean0, null())" --o

Limit to Mouginot 2019 mask

Actually, limit to approximate Mouginot 2019 mask - its a bit narrow in some places

# r.mask mask_GIC@Mouginot_2019 --o
r.grow input=mask_GIC@Mouginot_2019 output=mask_GIC_Mouginot_2019_grow radius=4.5 # three cells
r.mask mask_GIC_Mouginot_2019_grow --o
r.mapcalc "gates_xy_clean2 = gates_xy_clean1" --o
r.mask -r

# r.univar map=gates_xy_clean1
# r.univar map=gates_xy_clean2

Remove gates in areas from manually-drawn KML mask

See ./dat/remove_gates_manual.kml

v.import input=./dat/remove_gates_manual.kml output=remove_gates_manual --o
r.mask -i vector=remove_gates_manual --o
r.mapcalc "gates_xy_clean3 = gates_xy_clean2" --o
r.mask -r

r.univar map=gates_xy_clean2
r.univar map=gates_xy_clean3

Final Gates

g.copy "gates_xy_clean3,gates_final" --o

Add meta-data to gates

Add:

Gate ID
Calculate the average x,y of the gate, and then from that ONE point, determine the following. Do this from the average point rather than for each gate pixel because some gates span multiple sectors, or different ends of the gate are nearer different names, etc.
- Average lon,lat of gate
- Nearest citet:mouginot_2019_glacier region, sector, and name
- Nearest citet:bjork_2015_brief name

Do this for both the area vector and the point vector so that we can export

KML and GeoPackage with gates and metadata
simple CSV w/ gates and metadata.

Gate ID

# db.droptable -f table=gates_final
# db.droptable -f table=gates_final_pts

# areas (clusters of gate pixels, but diagonals are separate)
r.to.vect input=gates_final output=gates_final type=area --o
v.db.dropcolumn map=gates_final column=label
v.db.dropcolumn map=gates_final column=value
v.db.addcolumn map=gates_final columns="gate INT"
v.what.rast map=gates_final raster=gates_gateID column=gate type=centroid

# # points (each individual gate pixel)
# r.to.vect input=gates_final output=gates_final_pts type=point --o
# v.db.dropcolumn map=gates_final_pts column=label
# v.db.dropcolumn map=gates_final_pts column=value
# v.db.addcolumn map=gates_final_pts columns="gate INT"
# v.what.rast map=gates_final_pts raster=gates_gateID column=gate type=point

Mean x,y

# v.db.addcolumn map=gates_final columns="x DOUBLE PRECSION, y DOUBLE PRECISION, mean_x INT, mean_y INT, area INT"
v.db.addcolumn map=gates_final columns="mean_x INT, mean_y INT"
v.to.db map=gates_final option=coor columns=x,y units=meters
v.to.db map=gates_final option=area columns=area units=meters

for G in $(db.select -c sql="select gate from gates_final"|sort -n|uniq); do
  db.execute sql="UPDATE gates_final SET mean_x=(SELECT AVG(x) FROM gates_final WHERE gate == ${G}) where gate == ${G}"
  db.execute sql="UPDATE gates_final SET mean_y=(SELECT AVG(y) FROM gates_final WHERE gate == ${G}) where gate == ${G}"
done

v.out.ascii -c input=gates_final columns=gate,mean_x,mean_y | cut -d"|" -f4- | sort -n|uniq | v.in.ascii input=- output=gates_final_pts skip=1 cat=1 x=2 y=3 --o
v.db.addtable gates_final_pts
v.db.addcolumn map=gates_final_pts columns="gate INT"
v.db.update map=gates_final_pts column=gate query_column=cat

#v.db.addcolumn map=gates_final_pts columns="mean_x INT, mean_y INT"
v.to.db map=gates_final_pts option=coor columns=mean_x,mean_y units=meters

Here we have:

db.select table=gates_final|head -n10 # cat|gate|x|y|mean_x|mean_y
db.select table=gates_final_pts|head # cat|gate|mean_x|mean_y

Mean lon,lat

v.what.rast map=gates_final_pts raster=lon@PERMANENT column=lon
v.what.rast map=gates_final_pts raster=lat@PERMANENT column=lat

v.db.addcolumn map=gates_final columns="mean_lon DOUBLE PRECISION, mean_lat DOUBLE PRECISION"
for G in $(db.select -c sql="select gate from gates_final"|sort -n|uniq); do
    db.execute sql="UPDATE gates_final SET mean_lon=(SELECT lon FROM gates_final_pts WHERE gate = ${G}) where gate = ${G}"
    db.execute sql="UPDATE gates_final SET mean_lat=(SELECT lat FROM gates_final_pts WHERE gate = ${G}) where gate = ${G}"
done

Sector, Region, Names, etc.

Sector Number
Region Code
Nearest Sector or Glacier Name

v.db.addcolumn map=gates_final columns="sector INT"
v.db.addcolumn map=gates_final_pts columns="sector INT"
v.distance from=gates_final to=sectors@Mouginot_2019 upload=to_attr column=sector to_column=cat
v.distance from=gates_final_pts to=sectors@Mouginot_2019 upload=to_attr column=sector to_column=cat

v.db.addcolumn map=gates_final columns="region VARCHAR(2)"
v.db.addcolumn map=gates_final_pts columns="region VARCHAR(2)"
v.distance from=gates_final to=sectors@Mouginot_2019 upload=to_attr column=region to_column=SUBREGION1
v.distance from=gates_final_pts to=sectors@Mouginot_2019 upload=to_attr column=region to_column=SUBREGION1

v.db.addcolumn map=gates_final columns="Mouginot_2019 VARCHAR(99)"
v.db.addcolumn map=gates_final_pts columns="Mouginot_2019 VARCHAR(99)"
v.distance from=gates_final to=sectors@Mouginot_2019 upload=to_attr column=Mouginot_2019 to_column=NAME
v.distance from=gates_final_pts to=sectors@Mouginot_2019 upload=to_attr column=Mouginot_2019 to_column=NAME

v.db.addcolumn map=gates_final columns="Bjork_2015 VARCHAR(99)"
v.db.addcolumn map=gates_final_pts columns="Bjork_2015 VARCHAR(99)"
v.distance from=gates_final to=names@Bjork_2015 upload=to_attr column=Bjork_2015 to_column=name
v.distance from=gates_final_pts to=names@Bjork_2015 upload=to_attr column=Bjork_2015 to_column=name

v.db.addcolumn map=gates_final columns="Zwally_2012 INT"
v.db.addcolumn map=gates_final_pts columns="Zwally_2012 INT"
v.distance from=gates_final to=Zwally_2012@Zwally_2012 upload=to_attr column=Zwally_2012 to_column=cat_
v.distance from=gates_final_pts to=Zwally_2012@Zwally_2012 upload=to_attr column=Zwally_2012 to_column=cat_

v.db.addcolumn map=gates_final columns="Moon_2008 INT"
v.db.addcolumn map=gates_final_pts columns="Moon_2008 INT"
v.distance from=gates_final to=GlacierIDs@NSIDC_0642 upload=to_attr column=Moon_2008 to_column=GlacierID
v.distance from=gates_final_pts to=GlacierIDs@NSIDC_0642 upload=to_attr column=Moon_2008 to_column=GlacierID

v.db.addcolumn map=gates_final columns="Moon_2008_dist INT"
v.db.addcolumn map=gates_final_pts columns="Moon_2008_dist INT"
v.distance from=gates_final to=GlacierIDs@NSIDC_0642 upload=dist column=Moon_2008_dist
v.distance from=gates_final_pts to=GlacierIDs@NSIDC_0642 upload=dist column=Moon_2008_dist

v.db.addcolumn map=gates_final columns="n_pixels INT"
v.db.addcolumn map=gates_final_pts columns="n_pixels INT"
for G in $(db.select -c sql="select gate from gates_final"|sort -n|uniq); do
    db.execute sql="UPDATE gates_final SET n_pixels=(SELECT SUM(area)/(200*200) FROM gates_final WHERE gate = ${G}) where gate = ${G}"
    # now copy that to the average gate location (point) table
    db.execute sql="UPDATE gates_final_pts SET n_pixels = (SELECT n_pixels FROM gates_final WHERE gate = ${G}) WHERE gate = ${G}"
done

Clean up

db.dropcolumn -f table=gates_final column=area
# db.dropcolumn -f table=gates_final column=cat

Export as metadata CSV

mkdir -p out
db.select sql="SELECT gate,mean_x,mean_y,lon,lat,n_pixels,sector,region,Bjork_2015,Mouginot_2019,Zwally_2012,Moon_2008,Moon_2008_dist from gates_final_pts" separator=, | sort -n | uniq  > ./out/gate_meta.csv

Export Gates to KML

v.out.ogr input=gates_final output=./tmp/gates_final_${VELOCITY_CUTOFF}_${BUFFER_DIST}.kml format=KML --o
# open ./tmp/gates_final_${VELOCITY_CUTOFF}_${BUFFER_DIST}.kml

Sensitivity of results to gate distance and cutoff velocity

Run the gate detection algorithm at a variety of cutoff velocities and buffer distances

for VELOCITY_CUTOFF in 10 20 30 40 50 60 70 80 90 100 125 150; do
  for BUFFER_DIST in 1000 2000 3000 4000 5000 6000 7000 8000 9000; do
      . ./gate_IO.sh
  done
done

Heatmap

Compute

rm ./tmp/gate_test.dat
for M in $(g.mapset -l | tr ' ' '\n' | grep gates); do
  g.mapset ${M} --quiet
  r.mask -r --quiet
  g.region -d

  ### To generate this heatmap for a sub-region or single glacier, PRIOR to
  ### this in PERMANENT mapset, zoom in, then "set computational region extent from display",
  ### then save w/ "g.region save=SUBREGION --o", then in the loop zoom in on each mapset with:
  # g.region res=200 align=vel_baseline@MEaSUREs.0478 -pa
  # g.region region=SUBREGION@PERMANENT --o -pa

  r.mapcalc "vel_eff = (if(gates_x == 1, abs(vx_baseline@MEaSUREs.0478), 0) + if(gates_y == 1, abs(vy_baseline@MEaSUREs.0478), 0))" --o
  r.mask gates_final --o
  # frink "(m/yr * m) * m * kg/m^3 * 1E12 -> Gt/yr"

  ### simple D w/ "raw" unadjusted thickness
  r.mapcalc "tmp_D = (vel_eff * 200) * thickness@BedMachine * 917 / pow(10.0, 12)" --o

  ### D w/ adjusted thickness
  # in BedMachine mapset, from the OLS fits.summary()
  # r.mapcalc "thick_fit_adj=if(thickness<20,380*log(vel_baseline@MEaSUREs.0478,10)-380,thickness)"
  # r.mapcalc "tmp_D = (vel_eff * 200) * thick_fit_adj@BedMachine * 917 / pow(10.0, 12)" --o

  r.univar tmp_D | grep sum
  echo ${M} $(r.univar tmp_D | grep sum) >> ./tmp/gate_test.dat
  g.remove -f type=raster name=tmp_D
done

cat ./tmp/gate_test.dat | sort

Display

Effective Velocity

<<MSGS_pretty_print>>
<<GRASS_config>>

NOTDONE All Mapsets

Effective velocity (because gates may only be valid in x or y direction) was calculated for the sensitivity test above at one time, but for all gate locations (buffer dist and speed cutoff). Now we need to calculate it at all times

g.mapsets gates_vel_buf
r.mask -r

mapset=$(g.mapset -l | tr ' ' '\n' | grep -E gates_??_????| head -n1) # DEBUG
for mapset in $(g.mapset -l | tr ' ' '\n' | grep -E gates_??_????); do
  g.mapset ${mapset} --quiet

  g.region -d
  r.mapcalc "MASK = if((gates_x == 1) | (gates_y == 1), 1, null())" --o
  VX=$(g.list -m type=raster pattern=VX_????_??_?? mapset=* | head -n1) # DEBUG
  for VX in $(g.list -m type=raster pattern=VX_????_??_?? mapset=*); do
    VY=${VX/VX/VY}
    EX=${VX/VX/EX}
    EY=${VX/VX/EY}
    DATE=$(echo $VX | cut -d"_" -f2- | cut -d@ -f1)
    echo $DATE
    r.mapcalc "vel_eff_${DATE} = if(gates_x == 1, if(${VX} == -2*10^9, 0, abs(${VX})), 0) + if(gates_y == 1, if(${VY} == -2*10^9, 0, abs(${VY})), 0)"
    r.mapcalc "err_eff_${DATE} = if(gates_x == 1, if(${EX} == -2*10^9, 0, abs(${EX})), 0) + if(gates_y == 1, if(${EY} == -2*10^9, 0, abs(${EY})), 0)"
  done

  VX=$(g.list -m type=raster pattern=vx_????_??_?? mapset=Sentinel1 | head -n1) # DEBUG
  for VX in $(g.list -m type=raster pattern=vx_????_??_?? mapset=Sentinel1); do
    VY=${VX/vx/vy}
    EX=${VX/vx/ex}
    EY=${VX/vx/ey}
    DATE=$(echo $VX | cut -d"_" -f2- | cut -d@ -f1)
    echo $DATE
    r.mapcalc "vel_eff_${DATE} = 365 * (if(gates_x@gates_50_2500 == 1, if(isnull(${VX}), 0, abs(${VX}))) + if(gates_y == 1, if(isnull(${VY}), 0, abs(${VY}))))"
    r.mapcalc "err_eff_${DATE} = 365 * (if(gates_x@gates_50_2500 == 1, if(isnull(${EX}), 0, abs(${EX}))) + if(gates_y == 1, if(isnull(${EY}), 0, abs(${EY}))))"
  done
done

Just one velocity cutoff & buffer distance

g.mapsets -l

r.mask -r

MAPSET=gates_vel_buf

g.mapset MEaSUREs.0478
g.region -d
r.mapcalc "MASK = if((gates_x@${MAPSET} == 1) | (gates_y@${MAPSET} == 1), 1, null())" --o
dates=$(g.list type=raster pattern=VX_????_??_?? | cut -d"_" -f2-)
parallel --bar "r.mapcalc \"vel_eff_{1} = if(gates_x@${MAPSET} == 1, if(VX_{1} == -2*10^9, 0, abs(VX_{1})), 0) + if(gates_y@${MAPSET} == 1, if(VY_{1} == -2*10^9, 0, abs(VY_{1})), 0)\"" ::: ${dates}
parallel --bar "r.mapcalc \"err_eff_{1} = if(gates_x@${MAPSET} == 1, if(EX_{1} == -2*10^9, 0, abs(EX_{1})), 0) + if(gates_y@${MAPSET} == 1, if(EY_{1} == -2*10^9, 0, abs(EY_{1})), 0)\"" ::: ${dates}


g.mapset MEaSUREs.0481
g.region -d
r.mapcalc "MASK = if((gates_x@${MAPSET} == 1) | (gates_y@${MAPSET} == 1), 1, null())" --o
dates=$(g.list type=raster pattern=VX_????_??_?? | cut -d"_" -f2-)
parallel --bar "r.mapcalc \"vel_eff_{1} = if(gates_x@${MAPSET} == 1, if(isnull(VX_{1}), 0, abs(VX_{1})), 0) + if(gates_y@${MAPSET} == 1, if(isnull(VY_{1}), 0, abs(VY_{1})), 0)\"" ::: ${dates}
parallel --bar "r.mapcalc \"err_eff_{1} = if(gates_x@${MAPSET} == 1, if(isnull(EX_{1}), 0, abs(EX_{1})), 0) + if(gates_y@${MAPSET} == 1, if(isnull(EY_{1}), 0, abs(EY_{1})), 0)\"" ::: ${dates}


g.mapset MEaSUREs.0646
g.region -d
r.mapcalc "MASK = if((gates_x@${MAPSET} == 1) | (gates_y@${MAPSET} == 1), 1, null())" --o
dates=$(g.list type=raster pattern=VX_????_??_?? | cut -d"_" -f2-)
parallel --bar "r.mapcalc \"vel_eff_{1} = if(gates_x@${MAPSET} == 1, if(isnull(VX_{1}), 0, abs(VX_{1})), 0) + if(gates_y@${MAPSET} == 1, if(isnull(VY_{1}), 0, abs(VY_{1})), 0)\"" ::: ${dates}
parallel --bar "r.mapcalc \"err_eff_{1} = if(gates_x@${MAPSET} == 1, if(isnull(EX_{1}), 0, abs(EX_{1})), 0) + if(gates_y@${MAPSET} == 1, if(isnull(EY_{1}), 0, abs(EY_{1})), 0)\"" ::: ${dates}


g.mapset MEaSUREs.0731
g.region -d
r.mapcalc "MASK = if((gates_x@${MAPSET} == 1) | (gates_y@${MAPSET} == 1), 1, null())" --o
dates=$(g.list type=raster pattern=VX_????_??_?? | cut -d"_" -f2-)
parallel --bar "r.mapcalc \"vel_eff_{1} = if(gates_x@${MAPSET} == 1, if(isnull(VX_{1}), 0, abs(VX_{1})), 0) + if(gates_y@${MAPSET} == 1, if(isnull(VY_{1}), 0, abs(VY_{1})), 0)\"" ::: ${dates}
parallel --bar "r.mapcalc \"err_eff_{1} = if(gates_x@${MAPSET} == 1, if(isnull(EX_{1}), 0, abs(EX_{1})), 0) + if(gates_y@${MAPSET} == 1, if(isnull(EY_{1}), 0, abs(EY_{1})), 0)\"" ::: ${dates}

g.mapset Mouginot_pre2000
g.region -d
r.mapcalc "MASK = if((gates_x@${MAPSET} == 1) | (gates_y@${MAPSET} == 1), 1, null())" --o
VX=$(g.list type=raster pattern=vx_????_??_?? | head -n1) # DEBUG
for VX in $(g.list type=raster pattern=vx_????_??_??); do
  VY=${VX/vx/vy}
  DATE=$(echo $VX | cut -d"_" -f2-)
  echo $DATE
  r.mapcalc "vel_eff_${DATE} = if(gates_x@${MAPSET} == 1, if(isnull(${VX}), 0, abs(${VX})), 0) + if(gates_y@${MAPSET} == 1, if(isnull(${VY}), 0, abs(${VY})), 0)"
done

g.mapset MEaSUREs.0766
g.region -d
r.mapcalc "MASK = if((gates_x@${MAPSET} == 1) | (gates_y@${MAPSET} == 1), 1, null())" --o
dates=$(g.list type=raster pattern=VX_????_??_?? | cut -d"_" -f2-)
parallel --bar "r.mapcalc \"vel_eff_{1} = if(gates_x@${MAPSET} == 1, if(isnull(VX_{1}), 0, abs(VX_{1})), 0) + if(gates_y@${MAPSET} == 1, if(isnull(VY_{1}), 0, abs(VY_{1})), 0)\"" ::: ${dates}
parallel --bar "r.mapcalc \"err_eff_{1} = if(gates_x@${MAPSET} == 1, if(isnull(EX_{1}), 0, abs(EX_{1})), 0) + if(gates_y@${MAPSET} == 1, if(isnull(EY_{1}), 0, abs(EY_{1})), 0)\"" ::: ${dates}

g.mapset promice
g.region -d

dates=$(g.list type=raster pattern=vx_????_??_?? | cut -d"_" -f2-)

parallel --bar "r.mapcalc \"vel_eff_{1} = 365 * (if(gates_x@${MAPSET} == 1, if(isnull(vx_{1}), 0, abs(vx_{1})), 0) + if(gates_y@${MAPSET} == 1, if(isnull(vy_{1}), 0, abs(vy_{1})), 0))\"" ::: ${dates}

parallel --bar "r.mapcalc \"err_eff_{1} = 365 * (if(gates_x@${MAPSET} == 1, if(isnull(ex_{1}), 0, abs(ex_{1})), 0) + if(gates_y@${MAPSET} == 1, if(isnull(ey_{1}), 0, abs(ey_{1})), 0))\"" ::: ${dates}

# fix return code of this script so make continues
MSG_OK "vel_eff DONE"

Export all data to CSV

<<MSGS_pretty_print>>
<<GRASS_config>>

MSG_OK "Exporting..."
g.mapset PERMANENT
g.region -dp

MAPSET=gates_vel_buf

VEL_baseline="vel_baseline@MEaSUREs.0478 vx_baseline@MEaSUREs.0478 vy_baseline@MEaSUREs.0478 vel_err_baseline@MEaSUREs.0478 ex_baseline@MEaSUREs.0478 ey_baseline@MEaSUREs.0478"
VEL_0478=$(g.list -m mapset=MEaSUREs.0478 type=raster pattern=vel_eff_????_??_?? separator=space)
ERR_0478=$(g.list -m mapset=MEaSUREs.0478 type=raster pattern=err_eff_????_??_?? separator=space)
VEL_0481=$(g.list -m mapset=MEaSUREs.0481 type=raster pattern=vel_eff_????_??_?? separator=space)
ERR_0481=$(g.list -m mapset=MEaSUREs.0481 type=raster pattern=err_eff_????_??_?? separator=space)
VEL_0646=$(g.list -m mapset=MEaSUREs.0646 type=raster pattern=vel_eff_????_??_?? separator=space)
ERR_0646=$(g.list -m mapset=MEaSUREs.0646 type=raster pattern=err_eff_????_??_?? separator=space)
VEL_0731=$(g.list -m mapset=MEaSUREs.0731 type=raster pattern=vel_eff_????_??_?? separator=space)
ERR_0731=$(g.list -m mapset=MEaSUREs.0731 type=raster pattern=err_eff_????_??_?? separator=space)
VEL_0766=$(g.list -m mapset=MEaSUREs.0766 type=raster pattern=vel_eff_????_??_?? separator=space)
ERR_0766=$(g.list -m mapset=MEaSUREs.0766 type=raster pattern=err_eff_????_??_?? separator=space)
VEL_SENTINEL=$(g.list -m mapset=promice type=raster pattern=vel_eff_????_??_?? separator=space)
ERR_SENTINEL=$(g.list -m mapset=promice type=raster pattern=err_eff_????_??_?? separator=space)
VEL_MOUGINOT=$(g.list -m mapset=Mouginot_pre2000 type=raster pattern=vel_eff_????_??_?? separator=space)
THICK=$(g.list -m mapset=DEM type=raster pattern=DEM_???? separator=space)

LIST="lon lat err_2D gates_x@${MAPSET} gates_y@${MAPSET} gates_gateID@${MAPSET} sectors@Mouginot_2019 regions@Mouginot_2019 bed@BedMachine thickness@BedMachine surface@BedMachine ${THICK} ${VEL_baseline} ${VEL_0478}
${VEL_0481} ${VEL_0646} ${VEL_0731} ${VEL_0766} ${VEL_SENTINEL} ${VEL_MOUGINOT} errbed@BedMachine ${ERR_0478} ${ERR_0481} ${ERR_0646} ${ERR_0731} ${ERR_0766} ${ERR_SENTINEL}"

mkdir tmp/dat
r.mapcalc "MASK = if(gates_final@${MAPSET}) | if(mask_GIC@Mouginot_2019) | if(vel_err_baseline@MEaSUREs.0478) | if(DEM_2020@DEM)" --o
parallel --bar "if [[ ! -e ./tmp/dat/{1}.bsv ]]; then (echo x\|y\|{1}; r.out.xyz input={1}) > ./tmp/dat/{1}.bsv; fi" ::: ${LIST}
r.mask -r

# combine individual files to one mega csv
cat ./tmp/dat/lat.bsv | cut -d"|" -f1,2 | datamash -t"|" transpose > ./tmp/dat_100_5000_t.bsv
for f in ./tmp/dat/*; do
  cat $f | cut -d"|" -f3 | datamash -t"|" transpose >> ./tmp/dat_100_5000_t.bsv
done
cat ./tmp/dat_100_5000_t.bsv |datamash -t"|" transpose | tr '|' ',' > ./tmp/dat_100_5000.csv
rm ./tmp/dat_100_5000_t.bsv

Compute Errors

Velocity v Thickness Errors

Is velocity uncertainty important relative to thickness uncertainty?
- ANS: No, proportional velocity uncertainty is an order of magnitude less than thickness uncertainty.

import pandas as pd
import numpy as np

df = pd.read_csv("./tmp/dat_100_5000.csv")

thick = df['thickness@BedMachine'].copy()
thick[thick < 50] = 50  # IS THIS REASONABLE? IMPORTANT?
thick[thick == 0] = 1
thick_err = np.abs(df['errbed@BedMachine'].values)


vel = df['vel_baseline@MEaSUREs.0478']
vel_err = df['vel_err_baseline@MEaSUREs.0478']

print("Relative thickness error:")
print((thick_err/thick*100).describe())
print("\n\nRelative velocity error:")
print((vel_err/vel*100).describe())

# import matplotlib.pyplot as plt
# plt.clf()
# ax = plt.hexbin((thick_err/thick), (vel_err/vel), mincnt=1, bins='log', xscale='log', yscale='log')
# plt.xlabel('Relative Thickness Error [Log_10]')
# plt.ylabel('Relative Velocity Error [Log_10]')

Do velocity errors scale w/ velocity?

ANS: Not really.

import pandas as pd
import numpy as np
import datetime

# df = pd.read_csv("./tmp/dat_100_5000.csv")
# vel = df['vel_baseline@MEaSUREs.0478']
# vel_err = df['vel_err_baseline@MEaSUREs.0478']


vel = pd.read_csv("./tmp/dat_100_5000.csv", usecols=(lambda c: ('vel_eff' in c)))
vel.rename(columns=lambda c: datetime.datetime(int(c[8:12]), int(c[13:15]), int(c[16:18])), inplace=True)
vel.replace(0, np.nan, inplace=True)
# vel = vel.loc[:,vel.columns.year < 2018] # drop 2018
vel.sort_index(axis='columns', inplace=True)

vel_err = pd.read_csv("./tmp/dat_100_5000.csv", usecols=(lambda c: ('err_eff' in c)))
vel_err.rename(columns=lambda c: datetime.datetime(int(c[8:12]), int(c[13:15]), int(c[16:18])), inplace=True)
# vel_err = vel_err.loc[:,vel_err.columns.year < 2018] # drop 2018
for c in vel.columns:
    if c not in vel_err.columns:
        vel_err[c] = 0 
vel_err.sort_index(axis='columns', inplace=True) 


import matplotlib.pyplot as plt
plt.clf()
ax = plt.hexbin(vel.values.flatten()+1,
                vel_err.values.flatten()+1,
                mincnt=1,
                bins='log',
                xscale='log', yscale='log',
                gridsize=100)
plt.xlabel("Velocity [m yr$^{-1}$]")
plt.ylabel("Velocity Error [m yr$^{-1}$]")

Results (GIS)

from uncertainties import unumpy
import pandas as pd
import numpy as np

df_out = pd.DataFrame(columns=['Val','Unit'])
df_out.index.name = 'Method'

df = pd.read_csv("./tmp/dat_100_5000.csv")

thick = df['thickness@BedMachine']
# thick[thick < 50] = 50  # IS THIS REASONABLE? IMPORTANT?
# vel = df['vel_baseline@MEaSUREs.0478']
vel = np.abs(df['vx_baseline@MEaSUREs.0478'])*df['gates_x@gates_vel_buf'] + np.abs(df['vy_baseline@MEaSUREs.0478'])*df['gates_y@gates_vel_buf']
D = 200  * thick * vel * 917 / 1E12
df_out.loc['D'] = [np.sum(D), 'Gt']

err_thick = np.abs(df['errbed@BedMachine'].values)
# err_vel = np.abs(df['vel_err_baseline@MEaSUREs.0478'])
err_vel = np.abs(df['ex_baseline@MEaSUREs.0478'])*df['gates_x@gates_vel_buf'] + np.abs(df['ey_baseline@MEaSUREs.0478'])*df['gates_y@gates_vel_buf']

e_th = 200 * err_thick * vel * 917 / 1E12
df_out.loc['Err (Thickness)'] = [np.sum(e_th), 'Gt']
df_out.loc['Err (Thickness %)'] = [np.sum(e_th)/np.sum(D)*100, '%']

e_vel = 200 * thick * err_vel * 917 / 1E12
df_out.loc['Err (Velocity)'] = [np.sum(e_vel), 'Gt']
df_out.loc['Err (Velocity %)'] = [np.sum(e_vel)/np.sum(D)*100, '%']
df_out.loc['Err (Combined)'] = [np.sum(e_vel+e_th), 'Gt']
df_out.loc['Err (%)'] = [np.sum(e_vel+e_th)/np.sum(D)*100, '%']
# The above assumes everything is systematic/independent/correlated

# If errors are all random and uncorrelated:
t = unumpy.uarray(thick, err_thick)
v = unumpy.uarray(vel, err_vel)
e = np.sum(200 * t * v * 917 / 1E12)
df_out.loc['Random Errors'] = [e, 'Gt']
df_out.loc['Random Errors (%)'] = [e.s/e.n*100, '%']

df_out

Results (Mouginot 2019 Sector)

from uncertainties import unumpy
import pandas as pd
import numpy as np

df = pd.read_csv("./tmp/dat_100_5000.csv")

err_sector = pd.DataFrame(columns=['D', 'E', 'E%'])
err_sector.index.name = 'Sector'

sectors = np.unique(df['sectors@Mouginot_2019'].values)
for s in sectors:
    sub_s = df[df['sectors@Mouginot_2019'] == s]
    thick = sub_s['thickness@BedMachine']
    # vel = sub_s['vel_baseline@MEaSUREs.0478']
    vel = np.abs(sub_s['vx_baseline@MEaSUREs.0478']*sub_s['gates_x@gates_vel_buf']) + np.abs(sub_s['vy_baseline@MEaSUREs.0478']*sub_s['gates_y@gates_vel_buf'])
    D = 200  * thick * vel * 917 / 1E12
    err_thick = np.abs(sub_s['errbed@BedMachine'].values)
    # err_thick[np.where(err_thick < 50)] = 50  # IS THIS REASONABLE? IMPORTANT?
    e_th = 200 * err_thick * vel * 917 / 1E12
    err_sector.loc[s] = [np.sum(D), np.sum(e_th), np.round(np.sum(e_th),10)/np.round(np.sum(D),10)*100]

err_sector.loc['GIS'] = np.sum(err_sector, axis=0)
err_sector.loc['GIS']['E%'] = err_sector.loc['GIS']['E']/err_sector.loc['GIS']['D']*100

err_sector.to_csv('./tmp/err_sector_mouginot.csv')

err_sector.rename(columns = {'D':'D [Gt]',
                             'E':'Error [Gt]',
                             'E%':'Error [%]'}, inplace=True)

err_sector

Results (Mouginot 2019 Region)

from uncertainties import unumpy
import pandas as pd
import numpy as np

df = pd.read_csv("./tmp/dat_100_5000.csv")

err_region = pd.DataFrame(columns=['D','E', 'E%'])
err_region.index.name = 'Region'

regions = np.unique(df['regions@Mouginot_2019'].values)
for r in regions:
   sub_r = df[df['regions@Mouginot_2019'] == r]
   thick = sub_r['thickness@BedMachine']
   vel = np.abs(sub_r['vx_baseline@MEaSUREs.0478']*sub_r['gates_x@gates_vel_buf']) + np.abs(sub_r['vy_baseline@MEaSUREs.0478']*sub_r['gates_y@gates_vel_buf'])
   D = 200  * thick * vel * 917 / 1E12
   err_thick = np.abs(sub_r['errbed@BedMachine'].values)
   # err_thick[np.where(err_thick < 50)] = 50  # IS THIS REASONABLE? IMPORTANT?
   e_th = 200 * err_thick * vel * 917 / 1E12
   err_region.loc[r] = [np.sum(D), np.sum(e_th), np.round(np.sum(e_th),10)/np.round(np.sum(D),10)*100]

err_region.loc['GIS'] = np.sum(err_region, axis=0)
err_region.loc['GIS']['E%'] = err_region.loc['GIS']['E']/err_region.loc['GIS']['D']*100

err_region.to_csv('./tmp/err_region_mouginot.csv')

err_region.rename(columns = {'D':'D [Gt]', 
                         'E':'Error [Gt]',
                         'E%':'Error [%]'}, inplace=True)

err_region

Results (Gate)

from uncertainties import unumpy
import pandas as pd
import numpy as np

df = pd.read_csv("./tmp/dat_100_5000.csv")

err_gate = pd.DataFrame(columns=['D','E', 'E%'])
err_gate.index.name = 'Gate'

gates = np.unique(df['gates_gateID@gates_vel_buf'].values)
for g in gates:
    sub = df[df['gates_gateID@gates_vel_buf'] == g]
    thick = sub['thickness@BedMachine']
    vel = np.abs(sub['vx_baseline@MEaSUREs.0478'])*sub['gates_x@gates_vel_buf'] + np.abs(sub['vy_baseline@MEaSUREs.0478'])*sub['gates_y@gates_vel_buf']
    D = 200  * thick * vel * 917 / 1E12
    err_thick = np.abs(sub['errbed@BedMachine'].values)
    # err_thick[np.where(err_thick < 50)] = 50  # IS THIS REASONABLE? IMPORTANT?
    e_th = 200 * err_thick * vel * 917 / 1E12
    err_gate.loc[g] = [np.sum(D), np.sum(e_th), np.sum(e_th)/np.sum(D)*100]

err_gate.loc['GIS'] = np.sum(err_gate, axis=0)
err_gate.loc['GIS']['E%'] = err_gate.loc['GIS']['E']/err_gate.loc['GIS']['D']*100

gate_meta = pd.read_csv("./out/gate_meta.csv")
err_gate['name'] = ''
for g in err_gate.index.values:
    if (g == 'GIS'): continue
    if (sum(gate_meta.gate == g) == 0): continue
    err_gate.loc[g,'name'] = gate_meta[gate_meta.gate == g].Mouginot_2019.values[0]

err_gate.to_csv('./tmp/err_gate.csv')
err_gate.rename(columns = {'D':'D [Gt]', 
                           'E':'Error [Gt]',
                           'E%':'Error [%]'}, inplace=True),

err_gate

Raw data to discharge product

Load data

What columns are in the file?
Don’t show all the “vel_eff_YYYY_MM_DD” and “err_eff_YYYY_MM_DD” columns.

head -n1 ./tmp/dat_100_5000.csv | tr ',' '\n' | grep -v "vel_eff_*" | grep -v "err_eff_*" | grep -v "dh_*" | sort | uniq | tr '\n' '\t'
echo "also: dh_YYYY@elev, vel_eff_YYYY_MM_DD@various, etc."

Adjust thickness

Adjust “bad” thickness

Here we perform a few thickness adjustments:

300: All ice <= 20 m thick is assumed bad and set to the minimum “good” thickness in a gate if good exists, or 300 m if it does not exist
400: All ice <= 50 m thick is set to 400 m thick
fit: All ice <= 20 m thick is fit to the log10(thickness) v. log10(velocity) relationship, even though it is not a good fit.

For testing, gate clumps 9 (all bad) and 546 (some bad)

Adjust thickness w thickness v. velocity fit.

Table of thickness adjustments

Baseline discharge values for various thickness adjustments

Here we calculate:

D_baseline_th_noadj: Discharge with no thickness adjustment
D_baseline_th_300: The baseline discharge
D_baseline_th_400: The discharge assuming the aggressive thickness adjustment
D_baseline_th_fit: The discharge assuming the fitted thickness adjustment
D_baseline: The baseline discharge - picked from our favorite of the above. TBD

Map of where thickness adjustments occur

Temporal changes in thickness

Each pixel has DEM at YYYY. Assume 08-01 (DOY ~213)
Assume DEM_2020 is baseline
Linear interpolate between annual values
Take derivative to get dh/dt

Discharge

And more importantly and long-term, we calculate the following time series discharge products, using our preferred method (fill w/ 300 m):

D: Discharge at gate scale
D_err: The discharge error at gate scale
D_fill: The fill percentage for each gate at each point in time
D_sector: Same, but at Mouginot 2019 sector scale
D_sector_err
D_sector_fill
D_region: Same, but at Mouginot 2019 region scale
D_region_err
D_region_fill
D_all: Same, but all GIS
D_all_err
D_all_fill

import pandas as pd
import numpy as np
filled_D = pd.DataFrame(index=['A','B'], columns=['t1','t3','t4'], data=[[8,9,7],[4,2,1]])
fill = filled_D/filled_D
fill.loc['B','t3'] = np.nan

no_filled_D = filled_D * fill
# filled_weighted_D = filled_D / filled_D.sum()
no_filled_weighted_D = no_filled_D / filled_D.sum()

r = ((filled_D*fill)/filled_D.sum()).sum()
r.round(2)

SAVE & RESTORE STATE

%store D
%store D_err
%store fill
%store D_gates
%store D_gates_err
%store D_gates_fill_weight
%store D_sectors
%store D_sectors_err
%store D_sectors_fill_weight
%store D_regions
%store D_regions_err
%store D_regions_fill_weight
%store D_all
%store D_all_err
%store D_all_fill_weight
%store meta

%store -r

D = D.T['2000':].T
D_err = D_err.T['2000':].T
fill = fill.T['2000':].T
D_gates = D_gates.T['2000':].T
D_gates_err = D_gates_err.T['2000':].T
D_gates_fill_weight = D_gates_fill_weight.T['2000':].T
D_sectors = D_sectors.T['2000':].T
D_sectors_err = D_sectors_err.T['2000':].T
D_sectors_fill_weight = D_sectors_fill_weight.T['2000':].T
D_regions = D_regions.T['2000':].T
D_regions_err = D_regions_err.T['2000':].T
D_regions_fill_weight = D_regions_fill_weight.T['2000':].T
D_all = D_all.T['2000':].T
D_all_err = D_all_err.T['2000':].T
D_all_fill_weight = D_all_fill_weight.T['2000':].T

Export Data

Crop time series

STARTDATE='1986'
D_all = D_all.T[STARTDATE:].T
D_all_err = D_all_err.T[STARTDATE:].T
D_all_fill_weight = D_all_fill_weight.T[STARTDATE:].T
D_gates = D_gates.T[STARTDATE:].T
D_gates_err = D_gates_err.T[STARTDATE:].T
D_gates_fill_weight = D_gates_fill_weight.T[STARTDATE:].T
D_sectors = D_sectors.T[STARTDATE:].T
D_sectors_err = D_sectors_err.T[STARTDATE:].T
D_sectors_fill_weight = D_sectors_fill_weight.T[STARTDATE:].T
D_regions = D_regions.T[STARTDATE:].T
D_regions_err = D_regions_err.T[STARTDATE:].T
D_regions_fill_weight = D_regions_fill_weight.T[STARTDATE:].T
D_all = D_all.T[STARTDATE:].T
D_all_err = D_all_err.T[STARTDATE:].T
D_all_fill_weight = D_all_fill_weight.T[STARTDATE:].T

README

cat << EOF > ./out/README.txt
<<README>>
EOF

cat << EOF | cut -c2- >> ./out/README.txt 
 * Version

 This version of this README generated from git commit: $(git describe --always)
EOF

README for "Greenland Ice Sheet solid ice discharge from 1986 through March 2020"
 
Reference paper: doi:10.5194/essd-12-1367-2020  https://doi.org/10.5194/essd-12-1367-2020
Data Citations: doi:10.22008/promice/data/ice_discharge

Source: https://github.com/GEUS-PROMICE/ice_discharge

* Usage instructions:

When using any of the following data, you are required to cite the paper and the data set.

* Data Descriptions

Data sets released as part of this work include:
+ Discharge data
+ Gates
+ Surface Elevation Change
+ Code

Each are described briefly below.

** Discharge Data

This data set is made up of the following files

| Filename            | Description                                                            |
|---------------------+------------------------------------------------------------------------|
| GIS_D.csv           | Greenland Ice Sheet cumulative discharge by timestamp                  |
| GIS_err.csv         | Errors for GIS_D.csv                                                   |
| GIS_coverage.csv    | Coverage for GIS_D.csv                                                 |
| GIS.nc              | Discharge, errors, and covarge for GIS                                 |
| region_D.csv        | Regional discharge                                                     |
| region_err.csv      | Errors for region_D.csv                                                |
| region_coverage.csv | Coverage for region_D.csv                                              |
| region.nc           | Discharge, errors, and covarge for GIS regions                         |
| sector_D.csv        | Sector discharge                                                       |
| sector_err.csv      | Errors for sector_D.csv                                                |
| sector_coverage.csv | Coverage for sector_D.csv                                              |
| sector.nc           | Discharge, errors, and covarge for GIS sectors                         |
| gate_D.csv          | Gate discharge                                                         |
| gate_err.csv        | Errors for gate_D.csv                                                  |
| gate_coverage.csv   | Coverage for gate_D.csv                                                |
| gate.nc             | Discharge, errors, and covarge for GIS gates - including gate metadata |
|---------------------+------------------------------------------------------------------------|
| gate_meta.csv       | Metadata for each gate                                                 |

D and err data have units [Gt yr-1].
Coverage is in range [0, 1]

** Surface elevation change

See doi:10.22008/promice/data/DTU/surface_elevation_change/v1.0.0

** Code

See https://github.com/GEUS-PROMICE/ice_discharge

Gates

D, err, coverage

D_gatesT = D_gates.T
D_gates_errT = D_gates_err.T
D_gates_fill_weightT = D_gates_fill_weight.T

D_gatesT.index.name = "Date"
D_gates_errT.index.name = "Date"
D_gates_fill_weightT.index.name = "Date"

D_gatesT.to_csv('./out/gate_D.csv')
D_gates_errT.to_csv('./out/gate_err.csv')
D_gates_fill_weightT.to_csv('./out/gate_coverage.csv')

Sectors

# meta_sector = pd.DataFrame(index=meta.groupby('sectors').first().index)
# meta_sector['mean x'] = meta.groupby('sectors').mean()['x'].round().astype(int)
# meta_sector['mean y'] = meta.groupby('sectors').mean()['y'].round().astype(int)
# meta_sector['n gates'] = meta.groupby('sectors').count()['gates'].round().astype(int)
# meta_sector['region'] = meta.groupby('sectors').first()['regions']

D_sectorsT = D_sectors.T
D_sectors_errT = D_sectors_err.T
D_sectors_fill_weightT = D_sectors_fill_weight.T

D_sectorsT.index.name = "Date"
D_sectors_errT.index.name = "Date"
D_sectors_fill_weightT.index.name = "Date"

# meta_sector.to_csv('./out/sector_meta.csv')
D_sectorsT.to_csv('./out/sector_D.csv')
D_sectors_errT.to_csv('./out/sector_err.csv')
D_sectors_fill_weightT.to_csv('./out/sector_coverage.csv')

# meta_sector.head(10)

Regions

# meta_region = pd.DataFrame(index=meta.groupby('regions').first().index)
# meta_region['n gates'] = meta.groupby('regions').count()['gates'].round().astype(int)

D_regionsT = D_regions.T
D_regions_errT = D_regions_err.T
D_regions_fill_weightT = D_regions_fill_weight.T
D_regionsT.index.name = "Date"
D_regions_errT.index.name = "Date"
D_regions_fill_weightT.index.name = "Date"

# meta_region.to_csv('./out/region_meta.csv')
D_regionsT.to_csv('./out/region_D.csv')
D_regions_errT.to_csv('./out/region_err.csv')
D_regions_fill_weightT.to_csv('./out/region_coverage.csv')

# meta_region.head(10)

GIS

D_all.index.name = "Date"
D_all_err.index.name = "Date"
D_all_fill_weight.index.name = "Date"

D_all.to_csv('./out/GIS_D.csv', float_format='%.3f', header=["Discharge [Gt yr-1]"])
D_all_err.to_csv('./out/GIS_err.csv', float_format='%.3f', header=["Discharge Error [Gt yr-1]"])
D_all_fill_weight.to_csv('./out/GIS_coverage.csv', float_format='%.3f', header=["Coverage [unit]"])

Gates

g.mapset gates_vel_buf

v.out.ogr input=gates_final output=./out/gates.kml format=KML --o
(cd out; zip gates.kmz gates.kml; rm gates.kml)
v.out.ogr input=gates_final output=./out/gates.gpkg format=GPKG --o
v.out.ogr input=gates_final output=./out/gates.geojson format=GeoJSON --o

Elevation change

Done manually. See DOI

Code

Make sure this Org file is tidy enough…

Distribute

(cd out; zip -e /media/kdm/promicedata/ice_discharge/gates/gates.zip gates*)
(cd out; zip -e /media/kdm/promicedata/ice_discharge/d/D.zip D*csv)
cp ./out/README.txt /media/kdm/promicedata/ice_discharge/

zip -e /media/kdm/promicedata/ice_discharge/code/mankoff_et_al.zip ice_discharge.org

cp ${DATADIR}/Khan_2016/dhdt_1995-2015_GrIS.txt /media/kdm/promicedata/ice_discharge/surface_elevation_change

CSV to NetCDF

import pandas as pd
import xarray as xr
import numpy as np
import subprocess
import datetime

GIS

csvfile = 'GIS'

df_D = pd.read_csv('./out/' + csvfile + '_D.csv', index_col=0, parse_dates=True)
df_err = pd.read_csv('./out/' + csvfile + '_err.csv', index_col=0, parse_dates=True)
df_coverage = pd.read_csv('./out/' + csvfile + '_coverage.csv', index_col=0, parse_dates=True)

# ds = df_D.to_xarray()

# ds = xr.Dataset({'time': df_D})
ds = xr.Dataset()

ds["time"] = ("time", df_D.index)
ds["time"].attrs["long_name"] = "time of measurement"
ds["time"].attrs["standard_name"] = "time"
ds["time"].attrs["axis"] = "T"
ds["time"].attrs["cf_role"] = "timeseries_id"

ds["discharge"] = ("time", df_D['Discharge [Gt yr-1]'])
ds["discharge"].attrs["long_name"] = "Discharge"
ds["discharge"].attrs["standard_name"] = "land_ice_mass_tranport_due_to_calving_and_ice_front_melting"
ds["discharge"].attrs["units"] = "Gt yr-1"
ds["discharge"].attrs["coordinates"] = "time"

ds["err"] = ("time", df_err['Discharge Error [Gt yr-1]'])
ds["err"].attrs["long_name"] = "Error"
ds["err"].attrs["standard_name"] = "Uncertainty"
ds["err"].attrs["units"] = "Gt yr-1"
ds["err"].attrs["coordinates"] = "time"

ds["coverage"] = ("time", df_coverage['Coverage [unit]'])
ds["coverage"].attrs["long_name"] = "Coverage"
ds["coverage"].attrs["standard_name"] = "Coverage"
# ds["coverage"].attrs["units"] = "-"
ds["coverage"].attrs["coordinates"] = "time"

ds.attrs["featureType"] = "timeSeries"
ds.attrs["title"] = "Greenland discharge"
ds.attrs["summary"] = "Greenland discharge"
ds.attrs["keywords"] = "Greenland; Ice Discharge; Calving; Submarine Melt"
# ds.attrs["Conventions"] = "CF-1.8"
ds.attrs["source"] = "git commit: " + subprocess.check_output(["git", "describe", "--always"]).strip().decode('UTF-8')
# ds.attrs["comment"] = "TODO"
# ds.attrs["acknowledgment"] = "TODO"
# ds.attrs["license"] = "TODO"
# ds.attrs["date_created"] = datetime.datetime.now().strftime("%Y-%m-%d")
ds.attrs["creator_name"] = "Ken Mankoff"
ds.attrs["creator_email"] = "kdm@geus.dk"
ds.attrs["creator_url"] = "http://kenmankoff.com"
ds.attrs["institution"] = "GEUS"
# ds.attrs["time_coverage_start"] = "TODO"
# ds.attrs["time_coverage_end"] = "TODO"
# ds.attrs["time_coverage_resolution"] = "TODO"
ds.attrs["references"] = "10.22008/promice/ice_discharge"
ds.attrs["product_version"] = 2.0

# NOTE: Compression here does not save space
# comp = dict(zlib=True, complevel=5)
# encoding = {var: comp for var in ds.data_vars} # all
# encoding = {var: comp for var in ['time','coverage']} # some

ds.to_netcdf('./out/GIS.nc', mode='w')#, encoding=encoding)

Region

csvfile = 'region'

df_D = pd.read_csv('./out/' + csvfile + '_D.csv', index_col=0, parse_dates=True)
df_err = pd.read_csv('./out/' + csvfile + '_err.csv', index_col=0, parse_dates=True)
df_coverage = pd.read_csv('./out/' + csvfile + '_coverage.csv', index_col=0, parse_dates=True)

ds = xr.Dataset()

ds["time"] = (("time"), df_D.index)
ds["time"].attrs["cf_role"] = "timeseries_id"
ds["time"].attrs["standard_name"] = "time"
# ds["time"].attrs["units"] = "day of year"
# ds["time"].attrs["calendar"] = "julian"
ds["time"].attrs["axis"] = "T"

ds["region"] = (("region"), df_D.columns)
ds["region"].attrs["long_name"] = "Region"
ds["region"].attrs["standard_name"] = "N/A"
ds["region"].attrs["comment"] = "Regions from Mouginot (2019)"

ds["discharge"] = (("region", "time"), df_D.T.values)
ds["discharge"].attrs["long_name"] = "Discharge"
ds["discharge"].attrs["standard_name"] = "land_ice_mass_tranport_due_to_calving_and_ice_front_melting"
ds["discharge"].attrs["units"] = "Gt yr-1"
ds["discharge"].attrs["coordinates"] = "time region"

ds["err"] = (("region", "time"), df_err.T.values)
ds["err"].attrs["long_name"] = "Error"
ds["err"].attrs["standard_name"] = "Uncertainty"
ds["err"].attrs["units"] = "Gt yr-1"
ds["err"].attrs["coordinates"] = "time region"

ds["coverage"] = (("region", "time"), df_coverage.T.values)
ds["coverage"].attrs["long_name"] = "Coverage"
ds["coverage"].attrs["standard_name"] = "Coverage"
# ds["coverage"].attrs["units"] = "-"
ds["coverage"].attrs["coordinates"] = "time region"

# ds["lat"] = (("station"), meta.loc['lat'].astype(np.float32))
# #ds["lat"].attrs["coordinates"] = "station"
# ds["lat"].attrs["long_name"] = "latitude"
# ds["lat"].attrs["standard_name"] = "latitude"
# ds["lat"].attrs["units"] = "degrees_north"
# ds["lat"].attrs["axis"] = "Y"

# ds["lon"] = (("station"), meta.loc['lon'].astype(np.float32))
# #ds["lon"].attrs["coordinates"] = "station"
# ds["lon"].attrs["long_name"] = "longitude"
# ds["lon"].attrs["standard_name"] = "longitude"
# ds["lon"].attrs["units"] = "degrees_east"
# ds["lon"].attrs["axis"] = "X"

# ds["alt"] = (("station"), meta.loc['elev'].astype(np.float32))
# ds["alt"].attrs["long_name"] = "height_above_mean_sea_level"
# ds["alt"].attrs["standard_name"] = "altitude"
# # ds["alt"].attrs["long_name"] = "height above mean sea level"
# # ds["alt"].attrs["standard_name"] = "height"
# ds["alt"].attrs["units"] = "m"
# ds["alt"].attrs["positive"] = "up"
# ds["alt"].attrs["axis"] = "Z"

ds.attrs["featureType"] = "timeSeries"
ds.attrs["title"] = "Greenland discharge"
ds.attrs["summary"] = "Greenland discharge per region"
ds.attrs["keywords"] = "Greenland; Ice Discharge; Calving; Submarine Melt"
# ds.attrs["Conventions"] = "CF-1.8"
ds.attrs["source"] = "git commit: " + subprocess.check_output(["git", "describe", "--always"]).strip().decode('UTF-8')
# ds.attrs["comment"] = "TODO"
# ds.attrs["acknowledgment"] = "TODO"
# ds.attrs["license"] = "TODO"
# ds.attrs["date_created"] = datetime.datetime.now().strftime("%Y-%m-%d")
ds.attrs["creator_name"] = "Ken Mankoff"
ds.attrs["creator_email"] = "kdm@geus.dk"
ds.attrs["creator_url"] = "http://kenmankoff.com"
ds.attrs["institution"] = "GEUS"
# ds.attrs["time_coverage_start"] = "TODO"
# ds.attrs["time_coverage_end"] = "TODO"
# ds.attrs["time_coverage_resolution"] = "TODO"
ds.attrs["references"] = "10.22008/promice/ice_discharge"
ds.attrs["product_version"] = 2.0

comp = dict(zlib=True, complevel=9)
encoding = {var: comp for var in ds.data_vars} # all

ds.to_netcdf('./out/region.nc', mode='w', encoding=encoding)

Sectors

csvfile = 'sector'

df_D = pd.read_csv('./out/' + csvfile + '_D.csv', index_col=0, parse_dates=True)
df_err = pd.read_csv('./out/' + csvfile + '_err.csv', index_col=0, parse_dates=True)
df_coverage = pd.read_csv('./out/' + csvfile + '_coverage.csv', index_col=0, parse_dates=True)

ds = xr.Dataset()

ds["time"] = (("time"), df_D.index)
ds["time"].attrs["cf_role"] = "timeseries_id"
ds["time"].attrs["standard_name"] = "time"
# ds["time"].attrs["units"] = "day of year"
# ds["time"].attrs["calendar"] = "julian"
ds["time"].attrs["axis"] = "T"

ds["sector"] = (("sector"), df_D.columns)
ds["sector"].attrs["long_name"] = "Sector"
ds["sector"].attrs["standard_name"] = "N/A"
ds["sector"].attrs["comment"] = "Sectors from Mouginot (2019)"

ds["discharge"] = (("sector", "time"), df_D.T.values)
ds["discharge"].attrs["long_name"] = "Discharge"
ds["discharge"].attrs["standard_name"] = "land_ice_mass_tranport_due_to_calving_and_ice_front_melting"
ds["discharge"].attrs["units"] = "Gt yr-1"
ds["discharge"].attrs["coordinates"] = "time sector"

ds["err"] = (("sector", "time"), df_err.T.values)
ds["err"].attrs["long_name"] = "Error"
ds["err"].attrs["standard_name"] = "Uncertainty"
ds["err"].attrs["units"] = "Gt yr-1"
ds["err"].attrs["coordinates"] = "time sector"

ds["coverage"] = (("sector", "time"), df_coverage.T.values)
ds["coverage"].attrs["long_name"] = "Coverage"
ds["coverage"].attrs["standard_name"] = "Coverage"
# ds["coverage"].attrs["units"] = "-"
ds["coverage"].attrs["coordinates"] = "time sector"

# ds["lat"] = (("station"), meta.loc['lat'].astype(np.float32))
# #ds["lat"].attrs["coordinates"] = "station"
# ds["lat"].attrs["long_name"] = "latitude"
# ds["lat"].attrs["standard_name"] = "latitude"
# ds["lat"].attrs["units"] = "degrees_north"
# ds["lat"].attrs["axis"] = "Y"

# ds["lon"] = (("station"), meta.loc['lon'].astype(np.float32))
# #ds["lon"].attrs["coordinates"] = "station"
# ds["lon"].attrs["long_name"] = "longitude"
# ds["lon"].attrs["standard_name"] = "longitude"
# ds["lon"].attrs["units"] = "degrees_east"
# ds["lon"].attrs["axis"] = "X"

# ds["alt"] = (("station"), meta.loc['elev'].astype(np.float32))
# ds["alt"].attrs["long_name"] = "height_above_mean_sea_level"
# ds["alt"].attrs["standard_name"] = "altitude"
# # ds["alt"].attrs["long_name"] = "height above mean sea level"
# # ds["alt"].attrs["standard_name"] = "height"
# ds["alt"].attrs["units"] = "m"
# ds["alt"].attrs["positive"] = "up"
# ds["alt"].attrs["axis"] = "Z"

ds.attrs["featureType"] = "timeSeries"
ds.attrs["title"] = "Greenland discharge"
ds.attrs["summary"] = "Greenland discharge per sector"
ds.attrs["keywords"] = "Greenland; Ice Discharge; Calving; Submarine Melt"
# ds.attrs["Conventions"] = "CF-1.8"
ds.attrs["source"] = "git commit: " + subprocess.check_output(["git", "describe", "--always"]).strip().decode('UTF-8')
# ds.attrs["comment"] = "TODO"
# ds.attrs["acknowledgment"] = "TODO"
# ds.attrs["license"] = "TODO"
# ds.attrs["date_created"] = datetime.datetime.now().strftime("%Y-%m-%d")
ds.attrs["creator_name"] = "Ken Mankoff"
ds.attrs["creator_email"] = "kdm@geus.dk"
ds.attrs["creator_url"] = "http://kenmankoff.com"
ds.attrs["institution"] = "GEUS"
# ds.attrs["time_coverage_start"] = "TODO"
# ds.attrs["time_coverage_end"] = "TODO"
# ds.attrs["time_coverage_resolution"] = "TODO"
ds.attrs["references"] = "10.22008/promice/ice_discharge"
ds.attrs["product_version"] = 2.0

comp = dict(zlib=True, complevel=9)
encoding = {var: comp for var in ds.data_vars} # all

ds.to_netcdf('./out/sector.nc', mode='w', encoding=encoding)

Gates

csvfile = 'gate'

df_D = pd.read_csv('./out/' + csvfile + '_D.csv', index_col=0, parse_dates=True)
df_err = pd.read_csv('./out/' + csvfile + '_err.csv', index_col=0, parse_dates=True)
df_coverage = pd.read_csv('./out/' + csvfile + '_coverage.csv', index_col=0, parse_dates=True)

meta = pd.read_csv("./out/gate_meta.csv")

ds = xr.Dataset()

ds["time"] = (("time"), df_D.index)
ds["time"].attrs["cf_role"] = "timeseries_id"
ds["time"].attrs["standard_name"] = "time"
# ds["time"].attrs["units"] = "day of year"
# ds["time"].attrs["calendar"] = "julian"
ds["time"].attrs["axis"] = "T"

ds["gate"] = (("gate"), df_D.columns.astype(np.int32))
ds["gate"].attrs["long_name"] = "Gate"
ds["gate"].attrs["standard_name"] = "N/A"

ds["discharge"] = (("gate", "time"), df_D.T.values.astype(np.float32))
ds["discharge"].attrs["long_name"] = "Discharge"
ds["discharge"].attrs["standard_name"] = "land_ice_mass_tranport_due_to_calving_and_ice_front_melting"
ds["discharge"].attrs["units"] = "Gt yr-1"
ds["discharge"].attrs["coordinates"] = "time gate"

ds["err"] = (("gate", "time"), df_err.T.values.astype(np.float32))
ds["err"].attrs["long_name"] = "Error"
ds["err"].attrs["standard_name"] = "Uncertainty"
ds["err"].attrs["units"] = "Gt yr-1"
ds["err"].attrs["coordinates"] = "time gate"

ds["coverage"] = (("gate", "time"), df_coverage.T.values.astype(np.float32))
ds["coverage"].attrs["long_name"] = "Coverage"
ds["coverage"].attrs["standard_name"] = "Coverage"
# ds["coverage"].attrs["units"] = "-"
ds["coverage"].attrs["coordinates"] = "time gate"

ds["mean_x"] = (("gate"), meta.mean_x.astype(np.int32))
ds["mean_x"].attrs["long_name"] = "Mean x coordinate of gate in EPSG:3413"
ds["mean_x"].attrs["standard_name"] = "Mean x"

ds["mean_y"] = (("gate"), meta.mean_y.astype(np.int32))
ds["mean_y"].attrs["long_name"] = "Mean y coordinate of gate in EPSG:3413"
ds["mean_y"].attrs["standard_name"] = "Mean y"

ds["mean_lon"] = (("gate"), meta.lon.astype(np.float32))
ds["mean_lon"].attrs["long_name"] = "Mean lon coordinate of gate"
ds["mean_lon"].attrs["standard_name"] = "Longitude"

ds["mean_lat"] = (("gate"), meta.lat.astype(np.float32))
ds["mean_lat"].attrs["long_name"] = "Mean lat coordinate of gate"
ds["mean_lat"].attrs["standard_name"] = "Latitude"

ds["sector"] = (("gate"), meta.sector.astype(np.int32))
ds["sector"].attrs["long_name"] = "Mouginot 2019 sector containing gate"

ds["region"] = (("gate"), meta.region)
ds["region"].attrs["long_name"] = "Mouginot 2019 region containing gate"

ds["Zwally_2012"] = (("gate"), meta.Zwally_2012)
ds["Zwally_2012"].attrs["long_name"] = "Zwally 2012 sector containing gate"

ds["name_Bjørk"] = (("gate"), meta.Bjork_2015)
ds["name_Bjørk"].attrs["long_name"] = "Nearest name from Bjørk (2015)"

ds["name_Mouginot"] = (("gate"), meta.Mouginot_2019)
ds["name_Mouginot"].attrs["long_name"] = "Nearest name from Mouginot (2019)"

ds["ID_Moon"] = (("gate"), meta.Moon_2008)
ds["ID_Moon"].attrs["long_name"] = "Moon 2008 glacier ID (a.k.a NSIDC 0642)"

ds["ID_Moon_dist"] = (("gate"), meta.Moon_2008_dist)
ds["ID_Moon_dist"].attrs["long_name"] = "Moon 2008 glacier ID (a.k.a NSIDC 0642) distance"

ds.attrs["featureType"] = "timeSeries"
ds.attrs["title"] = "Greenland discharge"
ds.attrs["summary"] = "Greenland discharge per gate"
ds.attrs["keywords"] = "Greenland; Ice Discharge; Calving; Submarine Melt"
# ds.attrs["Conventions"] = "CF-1.8"
ds.attrs["source"] = "git commit: " + subprocess.check_output(["git", "describe", "--always"]).strip().decode('UTF-8')
# ds.attrs["comment"] = "TODO"
# ds.attrs["acknowledgment"] = "TODO"
# ds.attrs["license"] = "TODO"
# ds.attrs["date_created"] = datetime.datetime.now().strftime("%Y-%m-%d")
ds.attrs["creator_name"] = "Ken Mankoff"
ds.attrs["creator_email"] = "kdm@geus.dk"
ds.attrs["creator_url"] = "http://kenmankoff.com"
ds.attrs["institution"] = "GEUS"
# ds.attrs["time_coverage_start"] = "TODO"
# ds.attrs["time_coverage_end"] = "TODO"
# ds.attrs["time_coverage_resolution"] = "TODO"
ds.attrs["references"] = "10.22008/promice/ice_discharge"
ds.attrs["product_version"] = 2.0

comp = dict(zlib=True, complevel=9)
encoding = {var: comp for var in ds.data_vars} # all

ds.to_netcdf('./out/gate.nc', mode='w', encoding=encoding)

Figures

Overview

Top contributors

import pandas as pd

# annual average
df = pd.read_csv('./out/sector_D.csv', index_col=0, parse_dates=True)\
       .resample('1D')\
       .mean()\
       .interpolate(method='time')\
       .resample('A')\
       .mean()\

top_few = df.iloc[-1].sort_values(ascending=False).head(8)
# names = {}
# for m in mouginot_names: names[m] = ""
# names["JAKOBSHAVN_ISBRAE"] = "Sermeq Kujalleq (Jakobshavn Isbræ)"
# names["KOGE_BUGT_C"] = "(Køge Bugt C)"
# names["ZACHARIAE_ISSTROM"] = "Zachariae Isstrøm"
# names["RINK_ISBRAE"] = 'Kangilliup Sermia (Rink Isbræ)',
# names["NIOGHALVFJERDSFJORDEN"] = "(Nioghalvfjerdsfjorden)"

meta = pd.read_csv("./out/gate_meta.csv")

out = pd.DataFrame(index=top_few.index, columns=["D","x","y"])
for name in out.index:
    out.loc[name]["D"] = top_few[name]
    out.loc[name]["x"] = meta[meta["Mouginot_2019"] == name]["mean_x"].values[0]
    out.loc[name]["y"] = meta[meta["Mouginot_2019"] == name]["mean_y"].values[0]

out.to_csv('./tmp/overview_D.csv')
out

Prep data

grass -c ./G/fig_overview
g.region -d
r.mask -r


###
### vel
###
r.mask mask_GIC@Mouginot_2019 --o
r.mapcalc "vel = vel_baseline@MEaSUREs.0478" --o
r.mask -r

r.mapcalc "vel2 = if(vel < 100, null(), vel)" --o
# r.mapcalc "vel2 = if(vel2 > max(vel2)*0.5, max(vel2)*0.5, vel2)" --o
# r.colors -g map=vel2 color=oranges

cat << EOF | r.colors map=vel2 rules=-
0 255:255:255
99 255:255:255
100 255:235:215
500 255:128:0
15000 255:128:0
nv 255:255:255
default 255:255:255
EOF

g.copy vector=sectors@Mouginot_2019,sectors --o
g.copy vector=regions@Mouginot_2019,regions --o
g.copy vector=gates_final@gates_vel_buf,gates --o


### get baseline discharge at each gate
v.in.ascii input=./tmp/overview_D.csv output=overview_D separator=, skip=1 x=2 y=3 --o


###
### Ocean Basemap
###
r.mask -i mask@BedMachine
g.region res=1000
r.relief input=bed@BedMachine output=ocean_relief zscale=30  --o
# r.colors map=bed@BedMachine color=srtm_plus
r.colors -ne map=bed@BedMachine color=water
r.shade shade=ocean_relief color=bed@BedMachine output=ocean_shade --o
g.region -d


r.mask -r 
g.region res=2500
r.relief input=bed@BedMachine output=land_relief zscale=30  --o
r.mapcalc "land_relief = land_relief" --o
# r.colors map=bed@BedMachine color=srtm_plus
# r.shade shade=land_relief color=bed@BedMachine output=land_shade --o
g.region -d




g.region -d
r.mask -r
g.copy mask_ice@BedMachine,ice --o
# r.mapcalc "ice = if(isnull(ice), 2, 1)" --o
cat << EOF | r.colors map=ice rules=-
# 2 0:0:0
1 255:255:255
EOF



r.mask -r
d.mon start=wx0
d.erase
d.rast land_relief
d.rast ocean_shade
# d.rast ocean_relief
# d.rast land_shade
d.rast ice
d.rast vel2


## now patch them all to there R,G,B channels so we can produce a single ps.map raster.

g.region res=1000
r.mapcalc --o <<EOF
ocean.r = if(isnull(ocean_shade), null(), r#ocean_shade)
ocean.g = if(isnull(ocean_shade), null(), g#ocean_shade)
ocean.b = if(isnull(ocean_shade), null(), b#ocean_shade)
land.r = if(isnull(land_relief), null(), r#land_relief)
land.g = if(isnull(land_relief), null(), g#land_relief)
land.b = if(isnull(land_relief), null(), b#land_relief)
ice.r = if(isnull(ice), null(), r#ice)
ice.g = if(isnull(ice), null(), g#ice)
ice.b = if(isnull(ice), null(), b#ice)
vel.r = if(isnull(vel2), null(), r#vel2)
vel.g = if(isnull(vel2), null(), g#vel2)
vel.b = if(isnull(vel2), null(), b#vel2)
EOF

r.patch -s input=vel.r,ice.r,ocean.r,land.r output=rr --o
r.patch -s input=vel.g,ice.g,ocean.g,land.g output=gg --o
r.patch -s input=vel.b,ice.b,ocean.b,land.b output=bb --o

r.colors map=rr,gg,bb color=grey

d.erase
d.rgb -n red=rr green=gg blue=bb


g.region res=1000
cat << EOF | ps.map input=- output=./figs/psmap/BASE.ps --o
border n
rgb rr gg bb
EOF
g.region -d
# convert ${CONVERTOPTS} -transparent white ./figs/psmap/vel.ps ./figs/psmap/vel.png



# gates at each of the top 7 dischargers
# for row in $(db.select -c table=overview_D | tr ' ' '_' | tac); do
g.region -d
r.mapcalc "vel3 = vel2" --o

# r.colors map=vel3 color=oranges
cat << EOF | r.colors map=vel3 rules=-
0 255:255:255
99 255:255:255
100 255:235:215
10000 255:128:0
20000 255:128:0
nv 255:255:255
default 255:255:255
EOF



mkdir -p ./figs/psmap
id=$(db.select -c sql='select cat from overview_D'|head -n1) # DEBUG
for id in $(db.select -c sql='select cat from overview_D'); do
  xy=$(db.select -c sql="select int_1,int_2 from overview_D where cat == '${id}'")
  x=$(echo ${xy}|cut -d"|" -f1)
  y=$(echo ${xy}|cut -d"|" -f2)
  g.region w=$x s=$y e=$(( $x + 1 )) n=$(( $y + 1))
  g.region w=w-15000 s=s-15000 e=e+15000 n=n+15000 res=200
  g.region save=${id} --o

  cat << EOF | ps.map input=- output=./figs/psmap/${id}.ps --o
border n
raster vel3
vareas gates
  color black
  end
EOF
  convert -trim ./figs/psmap/${id}.{ps,png}
done


# overview map w/ boxes
# db.select -c sql='select cat from overview_D'
# 1 through 8

g.region -d res=1000
cat << EOF | ps.map input=- output=./figs/psmap/BASE_boxes.ps --o
border n
rgb rr gg bb
region 1
  color black
  width 1
  end  
region 2
  color black
  width 1
  end  
region 3
  color black
  width 1
  end  
region 4
  color black
  width 1
  end  
region 5
  color black
  width 1
  end  
region 6
  color black
  width 1
  end  
region 7
  color black
  width 1
  end  
region 8
  color black
  width 1
  end  
EOF
g.region -d

ps.map overview

create named regions with
- d.mon start=wx0
- d.rast vel
- d.rast gates color=red
- zoom, then “set computational region extent from display”
- then “g.region save=<name>”

g.mapset fig_overview
g.list type=region mapset=.
# g.remove -f type=region pattern=*

border n

rgb rr gg bb

# raster vel2

ps.map inset

border n

raster vel3

vareas gates
  color black
  end

end

Run psmap

g.region -d res=5000 -pa
ps.map input=overview.psmap output=overview.ps --o; o overview.ps
# ps.map -e input=overview.psmap output=overview.eps --o; o overview.eps

region=$(g.list type=region mapset=.|head -n1) # DEBUG
for region in $(g.list type=region mapset=.); do
  g.region ${region}
  g.region res=100
  ps.map input=overview_inset.psmap output=${region}.ps --o; 
  # o ${region}.ps
  convert -trim -transparent white ${region}.ps ./figs/psmap/${region}.png
  rm ${region}.ps
  # o ./figs/psmap/${region}.png
done

Configure in Inkscape

Manually align and adjust in Inkscape, add text, boxes, etc. then…

inkscape -z ./figs/overview.svg -e ./figs/overview.png

Discharge Time Series - GIS

import matplotlib.gridspec as gridspec
import numpy as np
import pandas as pd
from adjust_spines import adjust_spines as adj
import matplotlib.pyplot as plt
import datetime as dt

plt.close(1)

fig = plt.figure(1, figsize=(9,5)) # w,h
fig.clf()
grid = plt.GridSpec(2, 1, height_ratios=[1,5], hspace=0.1) # h, w

ax_D = fig.add_subplot(grid[1,:])
ax_coverage = fig.add_subplot(grid[0,:], sharex=ax_D)

adj(ax_D, ['left','bottom'])
adj(ax_coverage, ['left'])
ax_coverage.minorticks_off()
ax_coverage.tick_params(length=0, which='both', axis='x')

D = pd.read_csv("./out/GIS_D.csv", index_col=0, parse_dates=True)
err = pd.read_csv("./out/GIS_err.csv", index_col=0, parse_dates=True)
coverage = pd.read_csv("./out/GIS_coverage.csv", index_col=0, parse_dates=True)

THRESH = coverage.values.flatten() < 0.5
D[THRESH] = np.nan
err[THRESH] = np.nan
coverage[THRESH] = np.nan

# Add t0 and t_end so that graph covers a nice time span
def pad_df(df):
    df = df.append(pd.DataFrame(index=np.array(['1986-01-01']).astype('datetime64[ns]')), sort=True)
    idx = str(df.index.year.max())+'-12-31'
    df = df.append(pd.DataFrame(index=np.array([idx]).astype('datetime64[ns]')), sort=True)
    df = df.sort_index()
    return df

D = pad_df(D)
err = pad_df(err)
coverage = pad_df(coverage)

MS=4
ax_D.errorbar(err.index, D.values, fmt='o', mfc='none', mec='k', ms=MS)
for i in np.arange(D.values.size):
    if np.isnan(D.values[i]): continue
    alpha = coverage.values[i][0]
    if alpha < 0: alpha = 0
    if alpha > 1: alpha = 1
    ax_D.errorbar(D.index[i], D.values[i],
                  yerr=err.values[i], ecolor='grey',
                  alpha=alpha,
                  marker='o', ms=MS, mfc='k', mec='k')

# Take annual average from daily interpolated rather than the existing samples.
D_day_year = D.resample('1D').mean().interpolate(method='time',limit_area='inside',axis=0).resample('A').mean()
err_day_year = err.resample('1D').mean().interpolate(method='time',limit_area='inside',axis=0).resample('A').mean()

# No annual average if few sample
num_obs = D.resample('Y').count().values
D_day_year[num_obs < 4] = np.nan
err_day_year[num_obs < 4] = np.nan

Z=99
D_day_year.plot(drawstyle='steps', linewidth=3, ax=ax_D, alpha=0.75, color='orange', zorder=Z)

ax_D.legend("", framealpha=0)
ax_D.set_xlabel('Time [Years]')
ax_D.set_ylabel('Discharge [Gt yr$^{-1}$]')

import matplotlib.dates as mdates

ax_D.xaxis.set_major_locator(mdates.YearLocator())
ax_D.minorticks_off()
ax_D.xaxis.set_tick_params(rotation=-90) #, ha="left", rotation_mode="anchor")
for tick in ax_D.xaxis.get_majorticklabels():
    tick.set_horizontalalignment("left")

ax_D.set_xlim(D.index[0], D.index[-1])

###
### Coverage
###

ax_coverage.scatter(coverage.index, coverage.values*100,
                    color='k',
                    marker='.',
                    alpha=0.25)
                  # linewidth=3,
ax_coverage.set_ylim(45,105)
ax_coverage.set_yticks([50,100])
ax_coverage.spines['left'].set_bounds(ax_coverage.get_ylim()[0],100)
ax_coverage.set_ylabel('Coverage [%]')
plt.setp(ax_coverage.get_xticklabels(), visible=False)

plt.savefig('./figs/discharge_ts.png', transparent=False, bbox_inches='tight', dpi=300)

disp = pd.DataFrame(index = D_day_year.index.year,
                    data = {'D' : D_day_year.values.flatten(), 
                            'err' : err_day_year.values.flatten()})
disp.index.name = 'Year'
disp

Discharge Time Series - Regions

import matplotlib.gridspec as gridspec
import numpy as np
import pandas as pd
from adjust_spines import adjust_spines as adj
import matplotlib.pyplot as plt
import datetime as dt

from cycler import cycler
plt.rcParams['axes.prop_cycle'] = cycler('color', ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', \
                                                   '#9467bd', '#8c564b', '#e377c2', '#bcbd22', '#17becf'])

plt.close(1)

fig = plt.figure(1, figsize=(9,7)) # w,h
fig.clf()
# fig.set_tight_layout(True)
grid = plt.GridSpec(2, 1, height_ratios=[1,6], hspace=0.1) # h, w

ax_D = fig.add_subplot(grid[1,:])
ax_coverage = fig.add_subplot(grid[0,:], sharex=ax_D)

from adjust_spines import adjust_spines as adj
adj(ax_D, ['left','bottom'])
adj(ax_coverage, ['left'])
ax_coverage.minorticks_off()
ax_coverage.tick_params(length=0, which='both', axis='x')


D = pd.read_csv("./out/region_D.csv", index_col=0, parse_dates=True)
err = pd.read_csv("./out/region_err.csv", index_col=0, parse_dates=True)
coverage = pd.read_csv("./out/region_coverage.csv", index_col=0, parse_dates=True)

THRESH = coverage < 0.5
D[THRESH] = np.nan
err[THRESH] = np.nan
coverage[THRESH] = np.nan

def pad_df(df):
    df = df.append(pd.DataFrame(index=np.array(['1986-01-01']).astype('datetime64[ns]')), sort=True)
    idx = str(df.index.year.max())+'-12-31'
    df = df.append(pd.DataFrame(index=np.array([idx]).astype('datetime64[ns]')), sort=True)
    df = df.sort_index()
    return df

D = pad_df(D)
err = pad_df(err)
coverage = pad_df(coverage)

### Take annual average from daily interpolated rather than the existing samples.
D_day_year = D.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()
err_day_year=err.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()

# No annual average if few sample
num_obs = D.resample('Y').count().values
D_day_year[num_obs<=3] = np.nan
err_day_year[num_obs<=3] = np.nan

MS=4
Z=99
for r in D.columns:
    e = ax_D.errorbar(D[r].index, D[r].values, fmt='o', mfc='none', ms=MS)
    C = e.lines[0].get_color()
    D_day_year[r].plot(drawstyle='steps', linewidth=2, ax=ax_D,
                       color=C,
                       # color='orange'
                       alpha=0.75, zorder=Z)
    for i in np.arange(D.index.size):
        if np.isnan(D.iloc[i][r]): continue
        alpha = coverage.iloc[i][r]
        if alpha < 0: alpha = 0
        if alpha > 1: alpha = 1
        ax_D.errorbar(D.iloc[i].name, D.iloc[i][r],
                      yerr=err.iloc[i][r], ecolor='gray',
                      marker='o', ms=MS,
                      # mfc='k', mec='k',
                      color=C,
                      mfc=C, mec=C,
                      alpha=alpha)

    tx = pd.Timestamp(str(D[r].dropna().index[-1].year) + '-01-01') + dt.timedelta(days=380)
    ty = D_day_year[r].dropna().iloc[-1]
    # if r in ['CE', 'SW']: ty=ty-4
    # if r == 'NE': ty=ty+4
    # if r == 'NO': ty=ty-2
    ax_D.text(tx, ty, r, verticalalignment='center', horizontalalignment='left')

    # if r in ['CE','NE','SE']:
    ax_coverage.scatter(coverage.index, coverage[r]*100,
                        marker='.',
                        alpha=0.25,
                        color=C)

ax_coverage.set_ylabel('Coverage [%]')
ax_coverage.set_ylim(45,105)
ax_coverage.set_yticks([50,100])
ax_coverage.spines['left'].set_bounds(ax_coverage.get_ylim()[0],100)

import matplotlib.dates as mdates
ax_D.xaxis.set_major_locator(mdates.YearLocator())

# plt.legend()
ax_D.legend("", framealpha=0)
ax_D.set_xlabel('Time [Years]')
ax_D.set_ylabel('Discharge [Gt yr$^{-1}$]')
ax_D.set_xlim(D.index[0], D.index[-1])
ax_D.set_xticklabels(D.index.year.unique())
# ax_D.set_yscale('log')

plt.setp(ax_coverage.get_xticklabels(), visible=False)

ax_D.xaxis.set_tick_params(rotation=-90)
for tick in ax_D.xaxis.get_majorticklabels():
    tick.set_horizontalalignment("left")

plt.savefig('./figs/discharge_ts_regions.png', transparent=False, bbox_inches='tight', dpi=300)
# plt.savefig('./figs/discharge_ts_regions.pdf', transparent=True, bbox_inches='tight', dpi=300)

Err_pct = (err_day_year.values/D_day_year.values*100).round().astype(int).astype(str)
Err_pct[Err_pct.astype(float)<0] = 'NaN'
tbl = (D_day_year.round().fillna(value=0).astype(int).astype(str) + ' ('+Err_pct+')')
tbl.index = tbl.index.year.astype(str)
tbl.columns = [_ + ' (Err %)' for _ in tbl.columns]
tbl

Region & Catchment Contributions

2000 to present mean

D = pd.read_csv("./out/GIS_D.csv", index_col=0, parse_dates=True)
err = pd.read_csv("./out/GIS_err.csv", index_col=0, parse_dates=True)

D = D[D.index.year >= 2000]\
    .resample('1D',axis='rows')\
    .mean()\
    .interpolate(method='time',limit_area='inside')\
    .mean()

err = err[err.index.year >= 2000]\
      .resample('1D',axis='rows')\
      .mean()\
      .interpolate(method='time',limit_area='inside')\
      .mean()

print(D.values)
print(err.values)

Discharge time series

D = pd.read_csv("./out/GIS_D.csv", index_col=0, parse_dates=True)
err = pd.read_csv("./out/GIS_err.csv", index_col=0, parse_dates=True)

D = D.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()
err = err.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()

df = pd.DataFrame(index=D.index)
df['D'] = D
df['err'] = err

print("Min: ", df.loc[df['D'].idxmin])
print("Max: ", df.loc[df['D'].idxmin])
df

Contribution from SE

D = pd.read_csv("./out/region_D.csv", index_col=0, parse_dates=True)
err = pd.read_csv("./out/region_err.csv", index_col=0, parse_dates=True)

D = D.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()
err = err.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()

print("D SE described:\n", D['SE'].describe())
print("\nD SE err percent:", (err['SE']/D['SE']*100).describe())

df = pd.DataFrame(index=D.index)
df['SE'] = D['SE']
df['SEerr'] = err['SE']
df

Contribution from SE relative to the total

D = pd.read_csv("./out/region_D.csv", index_col=0, parse_dates=True)
err = pd.read_csv("./out/region_err.csv", index_col=0, parse_dates=True)

D = D.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()
err = err.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()

D_SE = D['SE']
# D_rest = D.drop(columns='SE').sum(axis='columns')
D_rest = D.sum(axis='columns')

print((D_SE/D_rest*100).describe())
print("")
print(D_SE/D_rest*100)

Contribution from others relative to the total

D = pd.read_csv("./out/region_D.csv", index_col=0, parse_dates=True)
err = pd.read_csv("./out/region_err.csv", index_col=0, parse_dates=True)

D = D.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()
err = err.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()

print("\nNO+NE+NW:", (D[['NO','NE','NW']].sum(axis='columns')).describe()[['min','max']])

# print("\nERR NO+NE+NW:", (err[['NO','NE','NW']].sum(axis='columns')).describe()[['mean']])

print("\nNO+NE+NW%:", (D[['NO','NE','NW']].sum(axis='columns')/D.sum(axis='columns')*100).describe()[['mean']])

print("\nNW increase:", D['NW'])
      
# D_SE = D['SE']
# # D_rest = D.drop(columns='SE').sum(axis='columns')
# D_rest = D.sum(axis='columns')

# print((D_SE/D_rest*100).describe())
# print("")
# print(D_SE/D_rest*100)

Sermeq Kujalleq (Jakobshavn Isbræ)

D = pd.read_csv("./out/sector_D.csv", index_col=0, parse_dates=True)
err = pd.read_csv("./out/sector_err.csv", index_col=0, parse_dates=True)

D = D['JAKOBSHAVN_ISBRAE'].resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()

err = err['JAKOBSHAVN_ISBRAE'].resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()

print("Jakobshavn described:\n", D.describe())

df = pd.DataFrame(index=D.index)
df['Jako'] = D
df['err'] = err
df

Discharge Time Series - Sectors

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from cycler import cycler
plt.rcParams['axes.prop_cycle'] = cycler('color', ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', \
                                                   '#9467bd', '#8c564b', '#e377c2', '#bcbd22', '#17becf'])

plt.close(1)

fig = plt.figure(1, figsize=(9,5)) # w,h
fig.clf()
# fig.set_tight_layout(True)
grid = plt.GridSpec(2, 1, height_ratios=[1,10], hspace=0.1) # h, w


ax_D = fig.add_subplot(grid[1,:])
ax_coverage = fig.add_subplot(grid[0,:], sharex=ax_D)

from adjust_spines import adjust_spines as adj
adj(ax_D, ['left','bottom'])
adj(ax_coverage, ['left'])
ax_coverage.minorticks_off()
ax_coverage.tick_params(length=0, which='both', axis='x')


D = pd.read_csv("./out/sector_D.csv", index_col=0, parse_dates=True)
err = pd.read_csv("./out/sector_err.csv", index_col=0, parse_dates=True)
coverage = pd.read_csv("./out/sector_coverage.csv", index_col=0, parse_dates=True)

THRESH = coverage < 0.5
D[THRESH] = np.nan
err[THRESH] = np.nan
coverage[THRESH] = np.nan


def pad_df(df):
    df = df.append(pd.DataFrame(index=np.array(['1986-01-01']).astype('datetime64[ns]')), sort=True)
    idx = str(df.index.year.max())+'-12-31'
    df = df.append(pd.DataFrame(index=np.array([idx]).astype('datetime64[ns]')), sort=True)
    df = df.sort_index()
    return df

D = pad_df(D)
err = pad_df(err)
coverage = pad_df(coverage)

### Take annual average from daily interpolated rather than the existing samples.
D_day_year = D.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()
err_day_year=err.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()


# No annual average if few sample
num_obs = D.resample('Y').count()
D_day_year[num_obs<=3] = np.nan
err_day_year[num_obs<=3] = np.nan


# largest average for last year
D_sort = D.resample('Y', axis='rows')\
          .mean()\
          .iloc[-1]\
          .sort_values(ascending=False)

LABELS={}
# for k in D_sort.head(8).index: LABELS[k] = k
# Use the last       ^ glaciers

# Make legend pretty
LABELS['JAKOBSHAVN_ISBRAE'] = 'Sermeq Kujalleq (Jakobshavn Isbræ)'
LABELS['HELHEIMGLETSCHER'] = 'Helheim Gletsjer'
LABELS['KANGERLUSSUAQ'] = 'Kangerlussuaq Gletsjer'
LABELS['KOGE_BUGT_C'] = '(Køge Bugt C)'
LABELS['ZACHARIAE_ISSTROM'] = 'Zachariae Isstrøm'
LABELS['RINK_ISBRAE'] = 'Kangilliup Sermia (Rink Isbræ)'
LABELS['NIOGHALVFJERDSFJORDEN'] = '(Nioghalvfjerdsbrae)'
LABELS['PETERMANN_GLETSCHER'] ='Petermann Gletsjer'

SYMBOLS={}
SYMBOLS['JAKOBSHAVN_ISBRAE'] = 'o'
SYMBOLS['HELHEIMGLETSCHER'] = 's'
SYMBOLS['KANGERLUSSUAQ'] = 'v'
SYMBOLS['KOGE_BUGT_C'] = '^'
SYMBOLS['NIOGHALVFJERDSFJORDEN'] = 'v'
SYMBOLS['RINK_ISBRAE'] = 's'
SYMBOLS['ZACHARIAE_ISSTROM'] = 'o'
SYMBOLS['PETERMANN_GLETSCHER'] ='^'

MS=4
Z=99
for g in LABELS.keys(): # for each glacier

    # scatter of symbols
    e = ax_D.errorbar(D.loc[:,g].index,
                      D.loc[:,g].values,
                      label=LABELS[g],
                      fmt=SYMBOLS[g],
                      mfc='none',
                      ms=MS)

    # Annual average
    C = e.lines[0].get_color()
    D_day_year.loc[:,g].plot(drawstyle='steps', linewidth=2,
                             label='',
                             ax=ax_D,
                             alpha=0.75, color=C, zorder=Z)

    # Error bars, each one w/ custom opacity.
    # Also, fill symbols w/ same opacity.
    for i,idx in enumerate(D.loc[:,g].index):
        alpha = coverage.loc[:,g].values[i]
        if np.isnan(alpha): alpha = 0
        if alpha < 0: alpha = 0
        if alpha > 1: alpha = 1
        ax_D.errorbar(D.loc[:,g].index[i],
                      D.loc[:,g].values[i],
                      yerr=err.loc[:,g].values[i],
                      alpha=alpha,
                      label='',
                      marker=SYMBOLS[g],
                      ecolor='grey',
                      color = C,
                      mfc=C, mec=C,
                      ms=MS)


    # Coverage scatter on top axis
    ax_coverage.scatter(D.loc[:,g].index,
                        coverage.loc[:,g].values*100,
                        alpha=0.25,
                        marker=SYMBOLS[g],
                        facecolor='none',
                        s=10,
                        color=C)

# yl = ax_D.get_ylim()

ax_D.legend(fontsize=8, ncol=2, loc=(0.0, 0.82), fancybox=False, frameon=False)
ax_D.set_xlabel('Time [Years]')
ax_D.set_ylabel('Discharge [Gt yr$^{-1}$]')
ax_D.set_xlim(D.index[0], D.index[-1])


import matplotlib.dates as mdates
ax_D.xaxis.set_major_locator(mdates.YearLocator())
ax_D.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
ax_D.xaxis.set_tick_params(rotation=-90)
for tick in ax_D.xaxis.get_majorticklabels():
    tick.set_horizontalalignment("left")

ax_coverage.set_ylabel('Coverage [%]')
ax_coverage.set_ylim([45,105])
ax_coverage.set_yticks([50,100])
ax_coverage.spines['left'].set_bounds(ax_coverage.get_ylim()[0],100)
plt.setp(ax_coverage.get_xticklabels(), visible=False)

plt.savefig('./figs/discharge_ts_topfew.png', transparent=False, bbox_inches='tight', dpi=300)

Err_pct = (err_day_year / D_day_year*100).round().fillna(value=0).astype(int).astype(str)
Err_pct = Err_pct[list(LABELS.keys())]
tbl = D_day_year[list(LABELS.keys())].round().fillna(value=0).astype(int).astype(str) + ' (' + Err_pct+')'
tbl.index = tbl.index.year.fillna(value=0).astype(str)
tbl.columns = [_ + ' (%)' for _ in tbl.columns]
tbl

Adjust Spines

http://matplotlib.org/examples/pylab_examples/spine_placement_demo.html

def adjust_spines(ax,spines, offset=10):
    for loc, spine in ax.spines.items():
        if loc in spines:
            spine.set_position(('outward', offset)) # outward by 10 points
            #spine.set_smart_bounds(True)
        else:
            spine.set_color('none') # don't draw spine

        # turn off ticks where there is no spine
        if 'left' in spines:
            ax.yaxis.set_tick_params(length=5)
            ax.yaxis.set_tick_params(direction='out')
            ax.yaxis.set_ticks_position('left')
            ax.yaxis.set_label_position('left')
        elif 'right' in spines:
            ax.yaxis.set_tick_params(length=5)
            ax.yaxis.set_tick_params(direction='out')
            ax.yaxis.set_ticks_position('right')
            ax.yaxis.set_label_position('right')
        else: # no yaxis ticks
            ax.yaxis.set_ticks([])

        if 'bottom' in spines:
            ax.xaxis.set_ticks_position('bottom')
            ax.xaxis.set_tick_params(length=5)
            ax.xaxis.set_tick_params(direction='out')
            ax.xaxis.set_label_position('bottom')
        elif 'top' in spines:
            ax.xaxis.set_ticks_position('top')
            ax.xaxis.set_tick_params(length=5)
            ax.xaxis.set_tick_params(direction='out')
            ax.xaxis.set_label_position('top')
        else: # no xaxis ticks
            ax.xaxis.set_ticks([])


if __name__ == "__main__":
    import numpy as np            
    x = np.random.random(100)
    fig = plt.figure(100)
    fig.clf()
    ax = fig.add_axes([0.1,0.1,0.8,0.8])
    ax.plot(x)
    adjust_spines(ax,["left","bottom"])

Auto update

Local

<<MSGS_pretty_print>>

MSG_OK "Updating Sentinel velocity files..."
workingdir=$(pwd)

# Warn the user that we need to manually update MEaSURES for now
MSG_WARN "MEaSUREs 0766.002 product needs to be manually updated if new data is available, see how in README.org in the datafolder"

# Update PROMICE IV 
cd ${DATADIR}/Promice200m/
if [[ -e urls.txt ]]; then cp urls.txt urls.txt.last; fi
curl https://dataverse.geus.dk/api/datasets/:persistentId?persistentId=doi:10.22008/FK2/LFZLZN | tr ',' '\n' | grep -E '"persistentId"' | cut -d'"' -f4 > urls.txt
chmod 777 urls.txt
if cmp -s urls.txt urls.txt.last; then
  MSG_WARN "No new Sentinel1 files..."
  #exit 255
fi

for PID in $(cat urls.txt | tail -n5); do
  wget --content-disposition --continue "https://dataverse.geus.dk/api/access/datafile/:persistentId?persistentId=${PID}"
done

MSG_OK "New Sentinel velocity files found..."
cd ${workingdir}

docker run --user $(id -u):$(id -g) --mount type=bind,src=${DATADIR},dst=/data --mount type=bind,src=$(pwd),dst=/home/user --env PARALLEL="--delay 0.1 -j -1" mankoff/ice_discharge:grass grass ./G/PERMANENT --exec ./update_worker.sh

cp ./tmp/dat_100_5000.csv ./tmp/dat_100_5000.csv.last

docker run --user $(id -u):$(id -g) --mount type=bind,src=${DATADIR},dst=/data --mount type=bind,src=$(pwd),dst=/home/user --env PARALLEL="--delay 0.1 -j -1" mankoff/ice_discharge:grass grass ./G/PERMANENT --exec ./export.sh

if cmp -s ./tmp/dat_100_5000.csv ./tmp/dat_100_5000.csv.last; then 
  MSG_ERR "No change"
  exit 255
fi

<<MSGS_pretty_print>>
<<GRASS_config>>
<<promice_import>>
<<MEaSURES_0766_import>>

MAPSET=gates_vel_buf
<<promice_effective_velocity>>
<<MEaSUREs_0766_effective_velocity>>

Upload

from pyDataverse.api import NativeApi
import os
import json

import subprocess
hash = subprocess.check_output(["git", "describe", "--always"]).strip().decode('utf-8')
hash = subprocess.check_output(["git", "describe", "--always", "--dirty='*'"]).strip().decode('utf-8')
assert("*" not in hash)

import secret
api_token = secret.dataverse_api_token

base_url = 'https://dataverse.geus.dk/'
api = NativeApi(base_url, api_token) # establish connection

# get dataverse metadata
identifier = 'doi:10.22008/promice/data/ice_discharge/d/v02'
resp = api.get_dataset(identifier)
files = resp.json()['data']['latestVersion']['files']
for f in files:
    persistentId = f['dataFile']['persistentId']
    description = f['dataFile']['description']
    filename = f['dataFile']['filename']
    fileId = f['dataFile']['id']
    print(filename)

    assert(os.path.isfile("./out/"+filename))

    description = description.split(".")[0] + ". "
    description = description + "Git hash: " + hash
    
    if 'content' in locals(): del(content)
    if filename[-3:] == ".nc": content = "application/x-netcdf"
    if filename[-3:] == "csv": content = "text/csv"
    if filename[-3:] == "txt": content = "text/plain"
    json_dict={"description":description, 
               "directoryLabel":".", 
               "forceReplace":True, 
               "filename":filename, 
               "label":filename, 
               "contentType":content}

    json_str = json.dumps(json_dict)
    d = api.replace_datafile(persistentId, "./out/"+filename, json_str)
  
    if d.json()["status"] == "ERROR": 
        print(d.content)
        print("\n")
        continue

    # need to update filenames after uploading because of DataVerse bug
    # https://github.com/IQSS/dataverse/issues/7223
    file_id = d.json()['data']['files'][0]['dataFile']['id']
    d2 = api.update_datafile_metadata(file_id, json_str=json_str, is_filepid=False)

resp = api.publish_dataset(identifier, "major")

Docker

GRASS

Dockerfile

FROM ubuntu:20.04

LABEL authors="Ken Mankoff"
LABEL maintainer="kdm@geus.dk"

# system environment
ENV DEBIAN_FRONTEND noninteractive

RUN apt-get -y update && apt-get install -y --no-install-recommends --no-install-suggests \
      bc \
      bsdmainutils \
      datamash \
      gdal-bin \
      grass \
      netcdf-bin \
      parallel \
      proj-bin \
      proj-data \
      zip \
  && apt-get autoremove -y \
  && apt-get clean -y \ 
  && rm -rf /var/lib/apt/lists/*

RUN echo LANG="en_US.UTF-8" > /etc/default/locale

ENV LANGUAGE en_US.UTF-8
ENV LANG C
ENV LC_ALL C
ENV LC_CTYPE C

ENV SHELL /bin/bash

# create a user
RUN useradd --create-home user && chmod a+rwx /home/user
ENV HOME /home/user
WORKDIR /home/user

RUN mkdir -p /data
ENV DATADIR /data

# switch the user
USER user

CMD ["/usr/bin/grass", "--version"]

Build

# docker build -f Dockerfile_grass -t ice_discharge_grass .
cd docker/grass
docker build -t ice_discharge_grass .
docker run -it ice_discharge_grass # run it

Test

container_args ?= run -it --cpus 7 --user $(shell id -u):$(shell id -g) --mount type=bind,src=$${DATADIR},dst=/data --mount type=bind,src=$(shell pwd),dst=/work --env PARALLEL="--delay 0.1" ice_discharge_grass

Deploy

# docker tag local-image:tagname new-repo:tagname
docker tag ice_discharge_grass mankoff/ice_discharge:grass
docker push mankoff/ice_discharge:grass

Python

Dockerfile and supporting files

FROM continuumio/miniconda3

RUN conda install \
  curl \
  cython \
  ipython \
  jupyter \
  matplotlib \
  numpy \
  pandas \
  pip \
  scipy \
  statsmodels \
  tabulate \
  xarray \
  && conda clean -a \
  && pip install --no-cache-dir \
  cfchecker \
  cfunits \
  grass-session \
  nc-time-axis \
  pyshp \
  semver \
  uncertainties \
  git+https://github.com/aussda/pyDataverse.git@3b040ff23b665ec2650bebcf4bd5478de6881af0

# create a user
RUN useradd --create-home user && chmod a+rwx /home/user
ENV HOME /home/user
WORKDIR /home/user

RUN mkdir -p /data
ENV DATADIR /data

# switch the user
USER user


# create a user
# RUN useradd -m -U user

# RUN chmod a+rwx /home/user
# ENV HOME /home/user
# RUN mkdir -p /data /work
# WORKDIR /work

# switch the user
# USER user

# RUN mkdir -p /data /work
# WORKDIR /work

# The code to run when container is started:
# ENTRYPOINT ["conda", "run", "-n", "ice_discharge", "python3"]

# For interactive shell
# RUN conda init bash
# RUN echo "conda activate ice_discharge" >> /root/.bashrc

Build

cd docker/conda
docker build -t ice_discharge_conda .
docker tag ice_discharge_conda:latest mankoff/ice_discharge:conda
docker run -it --mount type=bind,src=$(pwd),dst=/work mankoff/ice_discharge:conda python -c 'import pandas as pd; print(pd)'

Deploy

# docker tag local-image:tagname new-repo:tagname
docker tag ice_discharge_conda mankoff/ice_discharge:conda
docker push mankoff/ice_discharge:conda

enviroment.yml

It is more reproducible to use the Docker image mankoff/ice_discharge:conda, but for record-keeping sake, here is the Python environmnet.

docker run mankoff/ice_discharge:conda conda env export -n base

Emacs batch config

(require 'package)
(setq package-enable-at-startup nil)
(add-to-list 'package-archives '("melpa" . "http://melpa.org/packages/") t)
(add-to-list 'package-archives '("org" . "http://orgmode.org/elpa/") t)

(package-initialize)
;; Bootstrap `use-package'
(unless (package-installed-p 'use-package)
  (package-refresh-contents)
  (package-install 'use-package))
(setq package-enable-at-startup nil)

(defvar use-package-verbose t)
(require 'use-package)
(setq use-package-always-ensure t)
(setq use-package-always-defer nil)
(use-package diminish)
(require 'bind-key)   


(setenv "PATH" (concat "/home/ice/miniconda3/bin:" (getenv "PATH")))
(setq exec-path (split-string (getenv "PATH") ":"))

(use-package org
  :config
  (org-babel-do-load-languages
   'org-babel-load-languages
   '(;(emacs-lisp. t)   
     (shell . t)
     (calc . t)
     (python . t)
     (jupyter . t)
     ))

  (setq org-confirm-babel-evaluate nil) ;; don't ask to eval code
  (setq org-src-fontify-natively t)
  (setq org-export-use-babel t)
  (setq org-adapt-indentation nil)
  (setq org-edit-src-content-indentation 0)
  (setq org-src-preserve-indentation t))

(use-package jupyter
  :after org
  :init 
  ;;(add-hook 'org-mode-hook (lambda () (jupyter-org-interaction-mode)))
  :config 
  (setq jupyter-org-resource-directory "~/tmp/ob-jupyter-figs/")
  (add-hook 'org-babel-after-execute-hook 'org-display-inline-images 'append)
  :commands (jupyter-repl-history-previous jupyter-repl-history-next)
  :bind (:map jupyter-repl-mode-map
              ("C-n" . jupyter-repl-history-next)
              ("C-p" . jupyter-repl-history-previous))
  :custom-face
  (jupyter-repl-input-prompt ((t (:foreground "#000000"))))
  (jupyter-repl-output-prompt ((t (:foreground "#000000"))))
  (jupyter-repl-traceback ((t (:background "#FFFFFF"))))
)

Supplementary Material

Errors by gate sorted by total D, err, err %

Top by Discharge

import pandas as pd
import numpy as np

err_gate = pd.read_csv('./tmp/err_gate.csv', index_col=0)
err_gate.loc['GIS'] = np.nan
err_gate.dropna(inplace=True)
err_gate.rename(columns = {'D':'D [Gt]', 
                           'E':'Error [Gt]',
                           'E%':'Error [%]'}, inplace=True)
 
err_gate.sort_values('D [Gt]', inplace=True, ascending=False)
err_gate = err_gate.iloc[:25]
err_gate

Top by Error [Gt]

import pandas as pd
err_gate = pd.read_csv('./tmp/err_gate.csv', index_col=0)
err_gate.loc['GIS'] = np.nan
err_gate.dropna(inplace=True)
err_gate.rename(columns = {'D':'D [Gt]', 
                           'E':'Error [Gt]',
                           'E%':'Error [%]',
                           'Ew':'Error Weighted',
                           'Ew%':'Error Weighted %'}, inplace=True)

err_gate.sort_values('Error [Gt]', inplace=True, ascending=False)
err_gate = err_gate.iloc[:25]
err_gate

Top by Error [%]

import pandas as pd
err_gate = pd.read_csv('./tmp/err_gate.csv', index_col=0)
err_gate.loc['GIS'] = np.nan
err_gate.dropna(inplace=True)
err_gate.rename(columns = {'D':'D [Gt]', 
                           'E':'Error [Gt]',
                           'E%':'Error [%]',
                           'Ew':'Error Weighted',
                           'Ew%':'Error Weighted %'}, inplace=True)

err_gate.sort_values('Error [%]', inplace=True, ascending=False)
err_gate = err_gate.iloc[:25]
err_gate

Bottom by Error [%]

import pandas as pd
err_gate = pd.read_csv('./tmp/err_gate.csv', index_col=0)
err_gate.loc['GIS'] = np.nan
err_gate.dropna(inplace=True)
err_gate.rename(columns = {'D':'D [Gt]', 
                           'E':'Error [Gt]',
                           'E%':'Error [%]',
                           'Ew':'Error Weighted',
                           'Ew%':'Error Weighted %'}, inplace=True)

err_gate.sort_values('Error [%]', inplace=True, ascending=True)
err_gate = err_gate.iloc[:25]
err_gate

Annual averages from observations or linear interpolation

Compute the annual average for all GIS and some other glaciers using two methods, and compare the differences

Method 1: Annual average from observations
Method 2: Annual average after linearly interpolating observations to daily resolution

NOTE - do this with the “rawest” velocity product prior to any interpolation.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import datetime as dt

if 'th' not in locals():
    <<load_data>>
    <<adjust_thickness>>
    <<adjust_thickness_fit>>

# What is D where velocity is directly observed, no velocity filling?
D = (fill*vel).apply(lambda c: c * (th['fit'].values * 200), axis=0) * 917 / 1E12

ann = pd.DataFrame(index = D.resample("A", axis='columns').mean().columns)
ann['obs'] = D.sum(axis='rows').resample('A').mean().values
ann['daily'] = np.nan
for y in np.arange(1985,2018):
    ann.loc[ann.index.year == y, 'daily'] \
        = D\
        .T[D.columns.year == y]\
        .T\
        .sum(axis='rows')\
        .resample('1D')\
        .mean()\
        .interpolate(method="time")\
        .resample("A")\
        .mean()\
        .values
    
ann['diff'] = ann['obs'] - ann['daily']
ann['diff [%]'] = 100 - ann['daily']/ann['obs']*100

print(ann.describe())

ann

Køge Bugt Y2K

D = pd.read_csv("./out/sector_D.csv", index_col=0, parse_dates=True)
err = pd.read_csv("./out/sector_err.csv", index_col=0, parse_dates=True)

S = "KOGE_BUGT_C"
D = D[S].resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()

err = err[S].resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()

print("Køge Bugt C described:\n", D.describe())

df = pd.DataFrame(index=D.index)
df['KB'] = D
df['err'] = err
df

QA / QC / tests

Files used in this work

MEaSUREs

0478

File List

find ${DATADIR}/MEaSUREs/NSIDC-0478.002 | sed "s:^${DATADIR}::" | LC_ALL=C sort

File Metadata

somefile=$(find ${DATADIR}/MEaSUREs/NSIDC-0478.002 -name "*.tif"  | LC_ALL=C sort | grep ".*mosaic200.*vv.*" | head -n1)
md5sum ${somefile}
echo ""
gdalinfo -mm ${somefile}

md5sum hashes of all data

cd ${DATADIR}/MEaSUREs/NSIDC-0478.002
# awk 'NR % 5 == 0' input > output # print every 5th line
# head -n2 # top two files
find . -type f | LC_ALL=C sort | parallel --keep-order "md5sum {}"

0646

List of times, and count of folders at that time

find ${DATADIR}/MEaSUREs/NSIDC-0646.003/ -type f -name "*vx_v03.0.tif" | rev | cut -d"_" -f3 | rev | LC_ALL=C sort | uniq -c

0731

File List

Subset to vx

(cd ${DATADIR}/MEaSUREs/NSIDC-0731.005; find -name "*vx*.tif"| LC_ALL=C sort)

File Metadata

somefile=$(find ${DATADIR}/MEaSUREs/NSIDC-0731.005 -name "*vx*.tif" | LC_ALL=C sort | head -n1)
md5sum ${somefile}
echo ""
gdalinfo -mm ${somefile}

md5sum hashes of all data

cd ${DATADIR}/MEaSUREs/NSIDC-0731.005
# awk 'NR % 5 == 0' input > output # print every 5th line
# head -n2 # top two files
find . -type f -name "*vx*.tif" | LC_ALL=C sort | tail -n16| parallel --keep-order "md5sum {}"

Sentinel 1

File List

(cd ${DATADIR}/Promice200m; find . -name "*.nc" | LC_ALL=C sort)

File Metadata

somefile=$(find ${DATADIR}/Promice200m/ -name "*.nc" | LC_ALL=C sort | head -n1)
# md5sum ${somefile}
echo ""
ncdump -chs ${somefile}

md5sum Hashes

cd ${DATADIR}/Promice200m/
# awk 'NR % 5 == 0' input > output # print every 5th line
# head -n2 # top two files
find . -type f -name "*.nc" | LC_ALL=C sort | tail -n16 | parallel --keep-order "md5sum {}"

Results

GIS

FILE=./out/GIS_D.csv
cat ${FILE}

Regions

FILE=./out/region_D.csv
cat ${FILE}

Some arbitrary glaciers

FILE=./out/sector_D.csv

# GLLIST='Date|ACADEMY|SERMEQ_KUJ|BOWDOIN|BUGT|HELH|STORE|KANGER|ZACH|RINK|NIOGHALV|PETERMANN'
GLLIST='Date|JAKOBSH|HELH|STORE|ZACH' # |RINK'

# head -n1 ${FILE} \
#   | tr ',' '\n' \
#   | cat -n \
#   | grep -E ${GLLIST}

COLS=$(\
head -n1 <(cat ${FILE}) \
  | tr ',' '\n' \
  | cat -n \
  | grep -E ${GLLIST} \
  | cut -d$'\t' -f1 
)

cut -d, -f$(echo $COLS|tr ' ' ',') <(cat ${FILE})

LaTeX Setup

(add-to-list 'org-latex-classes
               `("copernicus"
                 "\\documentclass{copernicus}
               [NO-DEFAULT-PACKAGES]
               [NO-PACKAGES]
               [NO-EXTRA]"
                 ("\\section{%s}" . "\\section*{%s}")
                 ("\\subsection{%s}" . "\\subsection*{%s}")
                 ("\\subsubsection{%s}" . "\\subsubsection*{%s}")
                 ("\\paragraph{%s}" . "\\paragraph*{%s}")
                 ("\\subparagraph{%s}" . "\\subparagraph*{%s}"))
               )

;; (org-add-link-type
;;  "citet"  (lambda (key) (org-open-file cby-references-file t nil key))
;;  (lambda (path desc format)
;;    (cond
;;     ((eq format 'latex) (format "\\citet{%s}" path))
;;     ((eq format 'ascii) (format "%s" desc))
;;     )))
;; (org-add-link-type
;;  "citep"  (lambda (key) (org-open-file cby-references-file t nil key))
;;  (lambda (path desc format)
;;    (cond
;;     ((eq format 'latex) (format "\\citep{%s}" path))
;;     ((eq format 'ascii) (format "%s" desc))
;;     )))

(setq-local org-latex-title-command "")

bash	4.4.18-2ubuntu1.3
gdal-bin	3.0.4+dfsg-1~bionic0
grass-core	7.8.6-1~bionic1
nco	4.7.2-1
netcdf-bin	1:4.6.0-2build1
parallel	20161222-1
proj-bin	7.0.0-1~bionic0
sed	4.4-2

Files

ice_discharge.org

Latest commit

History

ice_discharge.org

File metadata and controls

Table of Contents

About This Document

Workflow

Summary

New in this version

Introduction

Input data

Methods

Terminology

Gate location

Thickness

Missing or invalid data

Invalid velocity

Missing velocity

Invalid thickness

Discharge

Discharge uncertainty

Results

Gates

Discharge

Discussion

Limitations

Data availability

Conclusions

Other

References

Figures

Overview

Heatmap

Ice Thickness v. Velocity 2D Histogram: Color = count

Discharge Time Series

Discharge Time Series: Regions

Discharge Time Series: Top Few

Appendix

Errors and uncertainties

Invalid thickness

Missing velocity

Errors from map projection

Velocity versus thickness

Software

Code

Misc Helper

Support pretty messages

GRASS config

Import Data

Bed and Surface

BedMachine v5

Sectors

Mouginot 2019

Import & Clean

Test

Zwally 2012

2D Area Error

Velocity

MEaSUREs

0478.002

Import

Baseline: Average of 2015-2017

Fill in holes

0481.004 TSX

Generate VRTs

Import VRTs

0646.003

Generate VRTs

Import VRTs

0731.005

Import

0766.002 Sentinel

Import

PROMICE IV 200m

Data Intro

Import data

Mouginot 2018 (pre-2000 velocities)

Display