/*----------------------------------------------------------------------*\ | Program to interpolate Peltier's topography grid on that of ETOPO5 | | | | This version uses bilinear interpolation because, and I quote, | | | | Bilinear interpolation is frequently "close enough for | | government work." -- Numerical Recipes in C, page 105 | | | | The interpolation calculates values from a one-degree grid whose top | | left cell is centered on (89.5N,-179.5W) to a five-minute grid whose | | top left cell is centered on (90N,180W). The geometry is shown by | | the diagram whose PostScript code is given below. In the diagram, | | solid lines delimit the one-degree grid cells, and open circles mark | | the centers of those cells. Dashed lines delimit the five-minute | | grid cells, and small filled dots mark the centers of those cells. | | The PostScript code for the diagram is as follows: | %! gsave 72 720 translate 200 dup scale 180 -90 translate 0 setlinewidth -180 90 moveto -180 88 lineto stroke -179 90 moveto -179 88 lineto stroke -178 90 moveto -178 88 lineto stroke -180 90 moveto -178 90 lineto stroke -180 89 moveto -178 89 lineto stroke -180 88 moveto -178 88 lineto stroke /s 1 12 div def /h s 2 div def [.01] 0 setdash /x1 -180 h sub def /x2 -178 h add def /y1 88 h sub def /y2 90 h add def x1 s x2 h add {dup y1 moveto y2 lineto stroke} for y1 s y2 h add {dup x1 exch moveto x2 exch lineto stroke} for [] 0 setdash -179.5 89.5 .02 0 360 arc stroke -179.5 88.5 .02 0 360 arc stroke -178.5 89.5 .02 0 360 arc stroke -178.5 88.5 .02 0 360 arc stroke /tinydot {0.005 0 360 arc fill} def /x1 -180 def /x2 -178 def /y1 88 def /y2 90 def y1 s y2 h add {/y exch def x1 s x2 h add {y tinydot} for } for showpage grestore % (end of PostScript code) | | | Peter N. Schweitzer (U.S. Geological Survey, Reston, VA 22092) | \*----------------------------------------------------------------------*/ #include #include #include #include #define WIDTH 360 #define HEIGHT 180 #define OFFSET 10000 #define DX 12 #define DY 12 static short *read_peltier_topography (char *filename) { int i,j,k,n; FILE *in; long size; short *d = NULL,*p; char *scratch; if (in = fopen (filename,"rb")) { if (fseek (in,0L,SEEK_END) == 0) { size = ftell (in); if (size != WIDTH * HEIGHT * sizeof(short)) { fprintf (stderr,"Error: expected %ld bytes, found %ld in %s\n", WIDTH*HEIGHT*sizeof(short),size,filename); exit (1); } rewind (in); } else size = WIDTH*HEIGHT*sizeof(short); if (p = d = (short *) malloc (size)) { if (scratch = (char *) malloc (WIDTH*sizeof(short))) { for (j=0; j < HEIGHT; j++) { if ((n = fread (scratch,sizeof(short),WIDTH,in)) == WIDTH) { swab (scratch,(char *)p,n*sizeof(short)); for (k=0; k < n; k++) *p++ -= OFFSET; } else { fprintf (stderr,"Error: only %ld of %ld words read from input file\n",n,WIDTH); exit (1); } } free (scratch); } else { fprintf (stderr,"Error: could not allocate space for scratch array\n"); exit (1); } } else { fprintf (stderr,"Error: could not allocate space for data from %s\n",filename); exit (1); } fclose (in); } else { fprintf (stderr,"Error: could not open input file %s\n",filename); exit (1); } return (d); } struct coefficients { double a; double b; double c; double d; }; static int calculate_coefficients (struct coefficients *z, int w, int h) { int i,j,k; double t; double u; k = 0; for (j=0; j < h; j++) { u = ((double)(h-j))/((double)h); for (i=0; i < w; i++) { t = ((double)i)/((double)w); z[k].a = (1.0 - t)*(1.0 - u); z[k].b = t*(1.0 - u); z[k].c = t*u; z[k].d = (1.0 - t)*u; k++; } } return (k); } main (int argc, char *argv[]) { char *input_file,*output_file; FILE *in,*out; long size; int i,j,k,m,n; struct coefficients z[DX*DY]; struct coefficients *pz; short *depth; short *p,*q; short *etopo5[DY]; short *e[DY]; double north_pole_depth; double south_pole_depth; long sum; double y1,y2,y3,y4; char *scratch; if (argc > 1) { input_file = argv[1]; depth = read_peltier_topography (input_file); if (argc > 2) { output_file = argv[2]; if (!(out = fopen (output_file,"wb"))) { fprintf (stderr,"Error: could not create output file %s\n",output_file); exit (1); } } else out = stdout; /*--------------------------------------------------------------*\ | Allocate DX rows of ETOPO5 grid cells and one scratch array. | \*--------------------------------------------------------------*/ size = DX*WIDTH*sizeof(short); for (j=0; j < DX; j++) if (etopo5[j] = (short *) malloc (size)) memset (etopo5[j],0,size); else { fprintf (stderr,"Error: could not allocate space for results\n"); exit (1); } if (!(scratch = (char *) malloc (size))) { fprintf (stderr,"Error: could not allocate scratch space\n"); exit (1); } /*--------------------------------------------------------------*\ | Calculate the depth at the north pole as the average of the | | top row of the Peltier data. | \*--------------------------------------------------------------*/ sum = 0L; p = depth; for (i=0; i < WIDTH; i++) sum += (long) *p++; north_pole_depth = ((double)sum)/(double)WIDTH; /*--------------------------------------------------------------*\ | Calculate the depth at the south pole as the average of the | | bottom row of the Peltier data. | \*--------------------------------------------------------------*/ sum = 0L; p = depth + WIDTH*(HEIGHT - 1); for (i=0; i < WIDTH; i++) sum += (long) *p++; south_pole_depth = ((double)sum)/(double)WIDTH; /*--------------------------------------------------------------*\ | Code for calculating the top six rows. | \*--------------------------------------------------------------*/ n = calculate_coefficients (z,DX,DY/2); q = depth; for (m=0; m < DY; m++) e[m] = etopo5[m]; /*----------------------------------------------------------*\ | Code for calculating the first six cells | \*----------------------------------------------------------*/ y1 = (double) q[WIDTH-1]; y2 = (double) q[0]; y3 = (double) north_pole_depth; y4 = (double) north_pole_depth; for (m=0; m < DY/2; m++) { for (k=DX/2; k < DX; k++) { pz = &z[m*DX + k]; *e[m]++ = (short) y1*pz->a + y2*pz->b + y3*pz->c + y4*pz->d; } } /*----------------------------------------------------------*\ | Code for calculating the middle cells | \*----------------------------------------------------------*/ for (i=0; i < WIDTH-1; i++) { y1 = (double) q[i]; y2 = (double) q[i+1]; y3 = (double) north_pole_depth; y4 = (double) north_pole_depth; for (m=0; m < DY/2; m++) { for (k=0; k < DX; k++) { pz = &z[m*DX + k]; *e[m]++ = (short) y1*pz->a + y2*pz->b + y3*pz->c + y4*pz->d; } } } /*----------------------------------------------------------*\ | Code for calculating the last six cells | \*----------------------------------------------------------*/ y1 = (double) q[WIDTH-1]; y2 = (double) q[0]; y3 = (double) north_pole_depth; y4 = (double) north_pole_depth; for (m=0; m < DY/2; m++) { for (k=0; k < DX/2; k++) { pz = &z[m*DX + k]; *e[m]++ = (short) y1*pz->a + y2*pz->b + y3*pz->c + y4*pz->d; } } /*----------------------------------------------------------*\ | Write the etopo5 arrays to the output file. | \*----------------------------------------------------------*/ for (m=0; m < DY/2; m++) { #ifdef NATIVE_ENDIAN_OUTPUT fwrite (etopo5[m],WIDTH*DX,sizeof(short),out); #else swab ((char *)etopo5[m],scratch,WIDTH*DX*sizeof(short)); fwrite (scratch,WIDTH*DX,sizeof(short),out); #endif } /*--------------------------------------------------------------*\ | Code for calculating the middle rows | \*--------------------------------------------------------------*/ n = calculate_coefficients (z,DX,DY); p = depth; for (j=0; j < HEIGHT-1; j++) { q = p + WIDTH; for (m=0; m < DY; m++) e[m] = etopo5[m]; /*----------------------------------------------------------*\ | Code for calculating the first six cells | \*----------------------------------------------------------*/ y1 = (double) q[WIDTH-1]; y2 = (double) q[0]; y3 = (double) p[0]; y4 = (double) p[WIDTH-1]; for (m=0; m < DY; m++) { for (k=DX/2; k < DX; k++) { pz = &z[m*DX + k]; *e[m]++ = (short) y1*pz->a + y2*pz->b + y3*pz->c + y4*pz->d; } } /*----------------------------------------------------------*\ | Code for calculating the middle cells | \*----------------------------------------------------------*/ for (i=0; i < WIDTH-1; i++) { y1 = (double) q[i]; y2 = (double) q[i+1]; y3 = (double) p[i+1]; y4 = (double) p[i]; for (m=0; m < DY; m++) { for (k=0; k < DX; k++) { pz = &z[m*DX + k]; *e[m]++ = (short) y1*pz->a + y2*pz->b + y3*pz->c + y4*pz->d; } } } /*----------------------------------------------------------*\ | Code for calculating the last six cells | \*----------------------------------------------------------*/ y1 = (double) q[i]; y2 = (double) q[0]; y3 = (double) p[0]; y4 = (double) p[i]; for (m=0; m < DY; m++) { for (k=0; k < DX/2; k++) { pz = &z[m*DX + k]; *e[m]++ = (short) y1*pz->a + y2*pz->b + y3*pz->c + y4*pz->d; } } /*----------------------------------------------------------*\ | Write the etopo5 arrays to the output file. | \*----------------------------------------------------------*/ for (m=0; m < DY; m++) { #ifdef NATIVE_ENDIAN_OUTPUT fwrite (etopo5[m],WIDTH*DX,sizeof(short),out); #else swab ((char *)etopo5[m],scratch,WIDTH*DX*sizeof(short)); fwrite (scratch,WIDTH*DX,sizeof(short),out); #endif } p += WIDTH; } /*--------------------------------------------------------------*\ | Code for calculating the bottom six rows | \*--------------------------------------------------------------*/ n = calculate_coefficients (z,DX,DY/2); for (m=0; m < DY; m++) e[m] = etopo5[m]; /*----------------------------------------------------------*\ | Code for calculating the first six cells | \*----------------------------------------------------------*/ y1 = (double) south_pole_depth; y2 = (double) south_pole_depth; y3 = (double) p[0]; y4 = (double) p[WIDTH-1]; for (m=0; m < DY/2; m++) { for (k=DX/2; k < DX; k++) { pz = &z[m*DX + k]; *e[m]++ = (short) y1*pz->a + y2*pz->b + y3*pz->c + y4*pz->d; } } /*----------------------------------------------------------*\ | Code for calculating the middle cells | \*----------------------------------------------------------*/ for (i=0; i < WIDTH-1; i++) { y1 = (double) south_pole_depth; y2 = (double) south_pole_depth; y3 = (double) p[i+1]; y4 = (double) p[i]; for (m=0; m < DY/2; m++) { for (k=0; k < DX; k++) { pz = &z[m*DX + k]; *e[m]++ = (short) y1*pz->a + y2*pz->b + y3*pz->c + y4*pz->d; } } } /*----------------------------------------------------------*\ | Code for calculating the last six cells | \*----------------------------------------------------------*/ y1 = (double) south_pole_depth; y2 = (double) south_pole_depth; y3 = (double) p[0]; y4 = (double) p[WIDTH-1]; for (m=0; m < DY/2; m++) { for (k=0; k < DX/2; k++) { pz = &z[m*DX + k]; *e[m]++ = (short) y1*pz->a + y2*pz->b + y3*pz->c + y4*pz->d; } } /*----------------------------------------------------------*\ | Write the etopo5 arrays to the output file. | \*----------------------------------------------------------*/ for (m=0; m < DY/2; m++) { #ifdef NATIVE_ENDIAN_OUTPUT fwrite (etopo5[m],WIDTH*DX,sizeof(short),out); #else swab ((char *)etopo5[m],scratch,WIDTH*DX*sizeof(short)); fwrite (scratch,WIDTH*DX,sizeof(short),out); #endif } } else { fprintf (stderr,"Usage: %s input_file output_file\n",argv[0]); exit (0); } return (0); } /*----------------------------------------------------------------------*\ \*----------------------------------------------------------------------*/