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.)

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-2016.
 * 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)
 *
 * Requires the 354MB etopo2.txt data file.
 *
 * Build: cc world.c -o world
 *
 */

#include 
#include 
#include 

#define FMT             "%22.4f "
#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 unsigned a,c,calif,d,h;
static int      doMountains,doOceans,mtnCriteria,seaCriteria,usFilter;
static double   all,depth,height,highest,landVolume,lowest,totalOceanV,v;
static double   arcticA,atlanticA,indianA,pacificA,southernA;
static double   arcticV,atlanticV,indianV,pacificV,southernV;
static FILE     *f;
static int      Atlantic(int a,int b);

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;
      }

  f = fopen("etopo2.txt","r");
  if (!f) {
    fprintf(stderr,"# Cannot find etopo2.txt data file.\n");
    fprintf(stderr,"# Here is sample output from a recent run:\n\n");
    fprintf(stderr,"         58320000.0000       data points\n");
    fprintf(stderr,"           989682.0000       coastal points (or 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,"         13235817.5995 km^2  Arctic Ocean area\n");
    fprintf(stderr,"         77366866.8857 km^2  Atlantic Ocean area\n");
    fprintf(stderr,"         63957480.5442 km^2  Indian Ocean area\n");
    fprintf(stderr,"        177939055.2334 km^2  Pacific Ocean area\n");
    fprintf(stderr,"         31387128.2118 km^2  Southern Ocean area\n");
    fprintf(stderr,"        363886348.4746 km^2  total ocean surface area\n");
    fprintf(stderr,"         17140286.9799 km^3  Arctic Ocean volume\n");
    fprintf(stderr,"        286869869.8852 km^3  Atlantic Ocean volume\n");
    fprintf(stderr,"        223806430.8616 km^3  Indian Ocean volume\n");
    fprintf(stderr,"        693164704.1121 km^3  Pacific Ocean volume\n");
    fprintf(stderr,"        110638551.7670 km^3  Southern Ocean volume\n");
    fprintf(stderr,"       1331619843.6058 km^3  total ocean volume\n");
    fprintf(stderr,"        117780843.0564 km^3  total land volume\n");
    fprintf(stderr,"       1331619843.6032 km^3  total ocean volume\n\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);
      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",(double)a);
  printf(FMT "      coastal points (or %.2f%%)\n",(double)c,(double)c/(double)d*100.0);
  printf(FMT "%%     coastal area\n",100.0*c/a);
  printf(FMT "%%     land area\n",100.0*h/a);
  printf(FMT "%%     water area\n",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",calif*3704/0.3048/5280);
  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 surface area\n",(arcticA+atlanticA+indianA+pacificA+southernA)/1e6);
  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);
  printf(FMT "km^3  total land volume\n",landVolume/1e9);
  printf(FMT "km^3  total ocean volume\n",totalOceanV/1e9);
  return 0;
}

int Atlantic(int a,int b)
{
  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: 19 Aug 2016