World Ocean Volumes

By Dan Allen

In 2004 I stumbled across some data that NASA posted. It was a dataset that for every two minutes of arc showed the average height of the earth or depth of the ocean. The data was called etopo and has since been revised a few times, and now there is a better dataset accurate to a single minute of arc. It is available here.

I realized that a variety of interesting things could be calculated from this data, so I wrote a program in C to parse the data. Here is the source code to my calculation too, but the 354 MB of data is too big for me to post. (Check the links given above.) Sometime I would like to adapt this to the new one minute of arc data.

Getting the byte order of binary data is really the only tricky thing about the program. In this case I had already written another program to print out the data in ASCII, so I can run the program everywhere without worrying about byte order.

The other tricky aspect to this program is deciding where the boundaries should be for the various oceans. This somewhat arbitrary aspect of such a program makes comparing results of other people hard to compare, because there are many legitimate ways to carve up the oceans. You can inspect my code below to see how I decided to do it. The Atlantic Ocean is the hardest to get right.

Embedded in the program are the results of running it from a previous run. That way I can see my results if I don't happen to have the 354 MB of data hanging around.

Enjoy!


/*
 * world.c - calculate volumes and areas of the world based on 2 minute squares.
 * world is Copyright Daniel K. Allen, 2004-2022.
 * All rights reserved.
 *
 *  2 Jul 2004 - Created by Dan Allen as etopo, 3:39:11 PM PDT.  (1,331,619,843,603,192,100 m^3)
 *  7 Jul 2004 - Volume of earth, average ht & depth. (Paradise,CA)
 * 14 Jul 2004 - Percent water, lat/lon at center of cells. (North Bend,WA)
 * 14 Dec 2006 - Renamed seas; calculates ocean areas and volumes. (Spring Lake,UT)
 * 15 Dec 2006 - Renamed world; added options; finally working for Atlantic ocean.
 * 19 Aug 2016 - Cleaned up. (Cape Cod,MA)
 * 13 Apr 2018 - Added int to main.
 * 26 Apr 2021 - Added land area, volume error; uses long double; runtime 5.73s on C50.
 * 10 Feb 2022 - Defaults to metric; also supports English units.
 *
 * Requires 353MB etopo2.txt data file from https://www.ngdc.noaa.gov/mgg/global/etopo2.html
 * Build: cc world.c -o world -Oz -mtune=skylake
 *
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#if defined(__APPLE__) || defined(applec) || defined(powerc)
  #define DEGSYM 161
#else
  #define DEGSYM 176
#endif

#define FMT             "%17.4Lf %-4s"
#define STARTLAT        (90*60-1)
#define STARTLON        (180*60)
#define D2R(x)          (x*3.14159265358979323846/180)
#define COS(x)          cos(D2R(x))
#define AREA(lat)       (3704*3704*fabs(COS((lat)/60.0)))
#define VOL(lat,depth)  (3704*3704*fabs(COS((lat)/60.0))*depth)

static int doMountains,doOceans,englishUnits,mtnCriteria,seaCriteria,usFilter;
static unsigned a,c,calif,d,h;
static long double all,depth,height,highest,landArea,landVolume,lowest,v;
static long double arcticA,atlanticA,indianA,pacificA,southernA;
static long double arcticV,atlanticV,indianV,pacificV,southernV,totalOceanV;
static long double M = 1,M2 = 1,M3 = 1;
static FILE *f;
static int Atlantic(int a,int b);

int main(int argc,char *argv[])
{
  int i,lat,lon,n;

  for (i = 1; i < argc; ++i)
    if (argv[i][0] == '-')
      switch (argv[i][1]) {
        case 'e': englishUnits = 1; M = 0.3048;
            M2 = 2.589988110336000293188; M3 = 4.168181825440579579833; break;
        case 'm': doMountains = 1; mtnCriteria = abs(atoi(argv[++i])); break;
        case 'o': doOceans = 1; seaCriteria = abs(atoi(argv[++i])); break;
        case 'u': usFilter = 1; break;
        default:
          fprintf(stderr,"# world [-e][-m heightMeters][-o depthMeters][-u]  # Calculate world areas & volumes.\n");
          return 1;
      }

  fprintf(stderr,"# 9 Jun 2004 etopo2.txt, 3888000 lines, 353808000 bytes\n");
  f = fopen("etopo2.txt","r");
  if (!f) {
    if (englishUnits) {
      fprintf(stderr,"    58320000.0000     data points\n");
      fprintf(stderr,"      989682.0000     coastal points (2.53%%)\n");
      fprintf(stderr,"           1.6970 %%   coastal area\n");
      fprintf(stderr,"          32.9255 %%   land area\n");
      fprintf(stderr,"          67.0745 %%   water area\n");
      fprintf(stderr,"        3759.0283 ft  average land height\n");
      fprintf(stderr,"       11089.9153 ft  average ocean depth\n");
      fprintf(stderr,"       -6200.8265 ft  average earth height & depth\n");
      fprintf(stderr,"       28192.2572 ft  highest land height\n");
      fprintf(stderr,"       34954.0682 ft  lowest ocean depth\n");
      fprintf(stderr,"        1077.1296 mi  California coastline\n");
      fprintf(stderr,"    56174376.4823 mi2 TOTAL land area\n");
      fprintf(stderr,"     5110377.7452 mi2 Arctic Ocean area\n");
      fprintf(stderr,"    29871524.3491 mi2 Atlantic Ocean area\n");
      fprintf(stderr,"    24694121.2930 mi2 Indian Ocean area\n");
      fprintf(stderr,"    68702643.2692 mi2 Pacific Ocean area\n");
      fprintf(stderr,"    12118637.9532 mi2 Southern Ocean area\n");
      fprintf(stderr,"   140497304.6097 mi2 TOTAL ocean area\n");
      fprintf(stderr,"    28257126.9654 mi3 TOTAL land volume\n");
      fprintf(stderr,"     4112173.5322 mi3 Arctic Ocean volume\n");
      fprintf(stderr,"    68823742.3167 mi3 Atlantic Ocean volume\n");
      fprintf(stderr,"    53694018.2157 mi3 Indian Ocean volume\n");
      fprintf(stderr,"   166299055.8337 mi3 Pacific Ocean volume\n");
      fprintf(stderr,"    26543600.1596 mi3 Southern Ocean volume\n");
      fprintf(stderr,"   319472590.0580 mi3 TOTAL ocean volume\n");
      fprintf(stderr,"        2.281e-07 mi3 TOTAL ocean volume error\n");
    }
    else {
      fprintf(stderr,"    58320000.0000     data points\n");
      fprintf(stderr,"      989682.0000     coastal points (2.53%%)\n");
      fprintf(stderr,"           1.6970 %%   coastal area\n");
      fprintf(stderr,"          32.9255 %%   land area\n");
      fprintf(stderr,"          67.0745 %%   water area\n");
      fprintf(stderr,"        1145.7518 m   average land height\n");
      fprintf(stderr,"        3380.2062 m   average ocean depth\n");
      fprintf(stderr,"       -1890.0119 m   average earth height & depth\n");
      fprintf(stderr,"        8593.0000 m   highest land height\n");
      fprintf(stderr,"       10654.0000 m   lowest ocean depth\n");
      fprintf(stderr,"        1733.4720 km  California coastline\n");
      fprintf(stderr,"   145490967.1946 km2 TOTAL land area\n");
      fprintf(stderr,"    13235817.5994 km2 Arctic Ocean area\n");
      fprintf(stderr,"    77366892.9017 km2 Atlantic Ocean area\n");
      fprintf(stderr,"    63957480.5441 km2 Indian Ocean area\n");
      fprintf(stderr,"   177939029.2160 km2 Pacific Ocean area\n");
      fprintf(stderr,"    31387128.2121 km2 Southern Ocean area\n");
      fprintf(stderr,"   363886348.4733 km2 TOTAL ocean area\n");
      fprintf(stderr,"   117780843.0562 km3 TOTAL land volume\n");
      fprintf(stderr,"    17140286.9799 km3 Arctic Ocean volume\n");
      fprintf(stderr,"   286869871.8832 km3 Atlantic Ocean volume\n");
      fprintf(stderr,"   223806430.8616 km3 Indian Ocean volume\n");
      fprintf(stderr,"   693164702.1142 km3 Pacific Ocean volume\n");
      fprintf(stderr,"   110638551.7672 km3 Southern Ocean volume\n");
      fprintf(stderr,"  1331619843.6060 km3 TOTAL ocean volume\n");
      fprintf(stderr,"        9.506e-07 km3 TOTAL ocean volume error\n");
    }
    return 1;
  }

  /* now read in 58,320,000 data points in 353,808,000 bytes */
  lat = STARTLAT;
  lon = STARTLON;
  while (fscanf(f,"%d",&n) == 1) {
    a++; all += n;
    if (n < 5 && n > -5) {
      c++;
      if (lat <= 41*60 && lat >= 32*60+32 && lon > 117*60 && lon < 125*60) calif++;
    }
    if (n > highest) highest = n;
    if (n < lowest) lowest = n;
    if (n < 0) {
      d++;
      depth += -n;
      v = VOL(lat,-n);
      totalOceanV += v;
      if (lat > 66*60) { /* N of the Arctic circle */
        arcticA += AREA(lat); arcticV += v;
      }
      else if (lat < -56*60) { /* S of Cape Horn */
        southernA += AREA(lat); southernV += v;
      }
      else if (lon < (-20-1)*60 && lon > (-113-1)*60) { /* Good Hope to West Australia */
        indianA += AREA(lat); indianV += v;
      }
      else if (Atlantic(lat,lon)) {
        atlanticA += AREA(lat); atlanticV += v;
      }
      else { /* West Australia to Cape Horn */
        pacificA += AREA(lat); pacificV += v;
      }

      if (doOceans && n <= -seaCriteria) { /* ocean depth in meters */
        if (usFilter && (lat < 27*60 || lat > 49*60 || lon < 60*60 || lon > 124*60)) goto skip;
        printf(FMT "%s%2d%c%02d'  %s%3d%c%02d'  %LG %s\n",
           v/1e9/M3,englishUnits?"mi3":"km3",
           lat>=0?"N":"S",abs(lat/60),DEGSYM,abs(lat%60),
           lon>=0?"W":"E",abs(lon-1)/60,DEGSYM,abs((lon-1)%60),
           n/M,englishUnits?"ft":"m");
      }
    }
    else {
      h++;
      height += n;
      v = VOL(lat,n);
      landArea += AREA(lat);
      landVolume += v;
      if (doMountains && n >= mtnCriteria) { /* mountain height in meters */
        if (usFilter && (lat < 27*60 || lat > 49*60 || lon < 60*60 || lon > 124*60)) goto skip;
        printf(FMT "%s%2d%c%02d'  %s%3d%c%02d'  %LG %s\n",
           v/1e9/M3,englishUnits?"mi3":"km3",
           lat>=0?"N":"S",abs(lat/60),DEGSYM,abs(lat%60),
           lon>=0?"W":"E",abs(lon-1)/60,DEGSYM,abs((lon-1)%60),
           n/M,englishUnits?"ft":"m");
      }
    }
    skip:
    lon -= 2;
    if (lon == -(STARTLON)) {
      lat -= 2;
      lon = STARTLON;
    }
  }
  fclose(f);

  /* print results */
  printf(FMT "data points\n",(long double)a,"");
  printf(FMT "coastal points (%.2Lf%%)\n",(long double)c,"",(long double)c/(long double)d*100.0);
  printf(FMT "coastal area\n",(long double)100.0*c/a,"%");
  printf(FMT "land area\n",(long double)100.0*h/a,"%");
  printf(FMT "water area\n",(long double)100.0*d/a,"%");
  printf(FMT "average land height\n",height/h/M,englishUnits?"ft":"m");
  printf(FMT "average ocean depth\n",depth/d/M,englishUnits?"ft":"m");
  printf(FMT "average earth height & depth\n",all/a/M,englishUnits?"ft":"m");
  printf(FMT "highest land height\n",highest/M,englishUnits?"ft":"m");
  printf(FMT "lowest ocean depth\n",-lowest/M,englishUnits?"ft":"m");
  printf(FMT "California coastline\n",(long double)calif*3704/M/(englishUnits?5280:1000),englishUnits?"mi":"km");
  printf(FMT "TOTAL land area\n",landArea/1e6/M2,englishUnits?"mi2":"km2");
  printf(FMT "Arctic Ocean area\n",arcticA/1e6/M2,englishUnits?"mi2":"km2");
  printf(FMT "Atlantic Ocean area\n",atlanticA/1e6/M2,englishUnits?"mi2":"km2");
  printf(FMT "Indian Ocean area\n",indianA/1e6/M2,englishUnits?"mi2":"km2");
  printf(FMT "Pacific Ocean area\n",pacificA/1e6/M2,englishUnits?"mi2":"km2");
  printf(FMT "Southern Ocean area\n",southernA/1e6/M2,englishUnits?"mi2":"km2");
  printf(FMT "TOTAL ocean area\n",(arcticA+atlanticA+indianA+pacificA+southernA)/1e6/M2,englishUnits?"mi2":"km2");
  printf(FMT "TOTAL land volume\n",landVolume/1e9/M3,englishUnits?"mi3":"km3");
  printf(FMT "Arctic Ocean volume\n",arcticV/1e9/M3,englishUnits?"mi3":"km3");
  printf(FMT "Atlantic Ocean volume\n",atlanticV/1e9/M3,englishUnits?"mi3":"km3");
  printf(FMT "Indian Ocean volume\n",indianV/1e9/M3,englishUnits?"mi3":"km3");
  printf(FMT "Pacific Ocean volume\n",pacificV/1e9/M3,englishUnits?"mi3":"km3");
  printf(FMT "Southern Ocean volume\n",southernV/1e9/M3,englishUnits?"mi3":"km3");
  printf(FMT "TOTAL ocean volume\n",(arcticV+atlanticV+indianV+pacificV+southernV)/1e9/M3,englishUnits?"mi3":"km3");
  v = fabsl(totalOceanV-(arcticV+atlanticV+indianV+pacificV+southernV))/1e9/M3;
  if (v) printf("%17.4Lg %-4sTOTAL ocean volume error\n",v,englishUnits?"mi3":"km3");
  return 0;
}

int Atlantic(int a,int b)
{
  long double lat,lon;
  lat = a / 60.0; lon = (b-1) / 60.0;
  /* between the Capes */
  if (lon < 69 && lon > -20) return 1;
  /* Mediterrean & Black Sea */
  if (lat >= 30 && lon < 0 && lon > -42) return 1;
  /* The Gulf of Mexico & Caribbean is tricky */
  if (lat > 0 && lon > 0 && lon < 100 && (lat*60 - 49.27)/(-60*lon) < -0.73) return 1;
  /* Otherwise not part of the Atlantic */
  return 0;
 }


 Created: 19 Aug 2016
Modified: 24 Jul 2024