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 accurage 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-2021.
 * 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.
 *
 * Requires 353MB etopo2.txt data file from https://www.ngdc.noaa.gov/mgg/global/etopo2.html
 * Build: cc world.c -o world -Os -mtune=skylake 
 *
 */

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

#define FMT             "%22.4Lf "
#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,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 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 '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 [-m height][-o depth][-u] # Calculate world areas & volumes.\n");
          fprintf(stderr,"# Note: West longitude positive, metric units used throughout.\n");
          return 1;
      }

  fprintf(stderr,"# 9 Jun 2004 etopo2.txt, 3888000 lines, 353808000 bytes\n");
  f = fopen("etopo2.txt","r");
  if (!f) {
    fprintf(stderr,"# Cached results from a recent run:\n\n");
    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,"        145490967.1946 km^2  total land area\n");
    fprintf(stderr,"         13235817.5994 km^2  Arctic Ocean area\n");
    fprintf(stderr,"         77366892.9017 km^2  Atlantic Ocean area\n");
    fprintf(stderr,"         63957480.5441 km^2  Indian Ocean area\n");
    fprintf(stderr,"        177939029.2160 km^2  Pacific Ocean area\n");
    fprintf(stderr,"         31387128.2121 km^2  Southern Ocean area\n");
    fprintf(stderr,"        363886348.4733 km^2  total ocean area\n");
    fprintf(stderr,"        117780843.0562 km^3  total land volume\n");
    fprintf(stderr,"         17140286.9799 km^3  Arctic Ocean volume\n");
    fprintf(stderr,"        286869871.8832 km^3  Atlantic Ocean volume\n");
    fprintf(stderr,"        223806430.8616 km^3  Indian Ocean volume\n");
    fprintf(stderr,"        693164702.1142 km^3  Pacific Ocean volume\n");
    fprintf(stderr,"        110638551.7672 km^3  Southern Ocean volume\n");
    fprintf(stderr,"       1331619843.6060 km^3  total ocean volume\n");
    fprintf(stderr,"           9.50625e-07 km^3  ocean volume error\n");
    return 1;
  }

  /* now read in 58,320,000 data points (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) { /* v > maxVol will print deep areas in oceans */
        if (usFilter && (lat < 27*60 || lat > 49*60 || lon < 60*60 || lon > 124*60)) goto skip;
        printf(FMT "m^3   %3dd %02d'  %4dd %02d'   %d m\n",
               v,lat/60,abs(lat%60),(lon-1)/60,abs((lon-1)%60),n);
      }
    }
    else {
      h++;
      height += n;
      v = VOL(lat,n);
      landArea += AREA(lat);
      landVolume += v;
      if (doMountains && n >= mtnCriteria) { /* v > 8e10 will print mountains */
        if (usFilter && (lat < 27*60 || lat > 49*60 || lon < 60*60 || lon > 124*60)) goto skip;
        printf(FMT "m^3   %3dd %02d'  %4dd %02d'   %d m\n",
               v,lat/60,abs(lat%60),(lon-1)/60,abs((lon-1)%60),n);
      }
    }
    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 (%.2f%%)\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 "ft    average land height\n",height/h/0.3048);
  printf(FMT "ft    average ocean depth\n",depth/d/0.3048);
  printf(FMT "ft    average earth height/depth\n",all/a/0.3048);
  printf(FMT "ft    highest land height\n",highest/0.3048);
  printf(FMT "ft    lowest ocean depth\n",-lowest/0.3048);
  printf(FMT "mi    California coastline\n",(long double)calif*3704/0.3048/5280);
  printf(FMT "km^2  total land area\n",landArea/1e6);
  printf(FMT "km^2  Arctic Ocean area\n",arcticA/1e6);
  printf(FMT "km^2  Atlantic Ocean area\n",atlanticA/1e6);
  printf(FMT "km^2  Indian Ocean area\n",indianA/1e6);
  printf(FMT "km^2  Pacific Ocean area\n",pacificA/1e6);
  printf(FMT "km^2  Southern Ocean area\n",southernA/1e6);
  printf(FMT "km^2  total ocean area\n",(arcticA+atlanticA+indianA+pacificA+southernA)/1e6);
  printf(FMT "km^3  total land volume\n",landVolume/1e9);
  printf(FMT "km^3  Arctic Ocean volume\n",arcticV/1e9);
  printf(FMT "km^3  Atlantic Ocean volume\n",atlanticV/1e9);
  printf(FMT "km^3  Indian Ocean volume\n",indianV/1e9);
  printf(FMT "km^3  Pacific Ocean volume\n",pacificV/1e9);
  printf(FMT "km^3  Southern Ocean volume\n",southernV/1e9);
  printf(FMT "km^3  total ocean volume\n",(arcticV+atlanticV+indianV+pacificV+southernV)/1e9);
  v = fabsl(totalOceanV-(arcticV+atlanticV+indianV+pacificV+southernV))/1e9;
  if (v) printf("%22Lg km^3  ocean volume error\n",v);
  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: 26 Apr 2021