merge standard release WRF/WPS V3.0.1.1 into wrffire
[wrffire.git] / WPS / geogrid / util / retile.c
blobdf055515ddeae2f04278340fe726f0e7c7b33a62
1 #include <stdio.h>
2 #include <sys/types.h>
3 #include <sys/stat.h>
4 #include <fcntl.h>
5 #include <unistd.h>
6 #include <math.h>
7 #include <string.h>
9 #define MSG_FLAG 0x80000001
11 #define N_BYTES 1
12 #define IN_TILE_DEGREES_LON 180.0
13 #define IN_TILE_DEGREES_LAT 180.0
14 #define IN_TILE_PTS_X 1250
15 #define IN_TILE_PTS_Y 1250
16 #define OUT_TILE_DEGREES_LON 180.0
17 #define OUT_TILE_DEGREES_LAT 180.0
18 #define OUT_TILE_PTS_X 1250
19 #define OUT_TILE_PTS_Y 1250
20 #define HALO_WIDTH 3
21 #define ZDIM 12
23 #define CACHESIZE 12
25 int **** data_cache = NULL;
26 char ** fname_cache = NULL;
27 int * lru = NULL;
29 int *** supertile = NULL;
30 int supertile_min_x = 999999;
31 int supertile_min_y = 999999;
33 void free_newtile(int *** x)
35 int i, j, k;
37 for(i=0; i<IN_TILE_PTS_X; i++)
39 for(j=0; j<IN_TILE_PTS_Y; j++)
41 free(x[i][j]);
43 free(x[i]);
45 free(x);
48 int *** read_tile(char * fname)
50 int *** retval;
51 int i, j, k, b;
52 unsigned char buf[IN_TILE_PTS_X*IN_TILE_PTS_Y*ZDIM*N_BYTES];
53 int fd;
55 retval = (int ***)malloc(sizeof(int **)*IN_TILE_PTS_X);
56 for(i=0; i<IN_TILE_PTS_X; i++)
58 retval[i] = (int **)malloc(sizeof(int *)*IN_TILE_PTS_Y);
59 for(j=0; j<IN_TILE_PTS_Y; j++)
61 retval[i][j] = (int *)malloc(sizeof(int)*ZDIM);
65 if ((fd = open(fname,O_RDONLY)) == -1)
67 fprintf(stderr,"Error opening source file %s\n", fname);
68 return 0;
71 read(fd, (void *)&buf, IN_TILE_PTS_X*IN_TILE_PTS_Y*ZDIM*N_BYTES);
73 close(fd);
75 /* Put buf into retval */
76 for(i=0; i<IN_TILE_PTS_X; i++)
78 for(j=0; j<IN_TILE_PTS_Y; j++)
80 for(k=0; k<ZDIM; k++)
82 retval[i][j][k] = 0;
83 for(b=0; b<N_BYTES; b++)
85 retval[i][j][k] |= buf[k*N_BYTES*IN_TILE_PTS_X*IN_TILE_PTS_Y+j*N_BYTES*IN_TILE_PTS_X+i*N_BYTES+b] << 8*(N_BYTES-b-1);
91 return retval;
94 int *** get_tile_from_cache(int i, int j)
96 int ii, jj, kk, k, least, least_idx;
97 int *** retval, *** localptr;
98 char * fname;
100 fname = (char *)malloc(256);
102 i = (i/IN_TILE_PTS_X)*IN_TILE_PTS_X+1;
103 j = (j/IN_TILE_PTS_Y)*IN_TILE_PTS_Y+1;
105 snprintf(fname,256,"%5.5i-%5.5i.%5.5i-%5.5i",i,i+IN_TILE_PTS_X-1,j,j+IN_TILE_PTS_Y-1);
107 /* Find out whether tile containing (i,j) is in cache */
108 if (data_cache != NULL)
110 for(k=0; k<CACHESIZE; k++)
112 if (fname_cache[k] != NULL)
114 if (strncmp(fname_cache[k],fname,256) == 0)
116 free(fname);
117 retval = data_cache[k];
118 return retval;
124 /* If not, read from file */
125 localptr = read_tile(fname);
127 /* Also store tile in the cache */
128 if (data_cache == NULL)
130 data_cache = (int ****)malloc(sizeof(int ***)*CACHESIZE);
131 fname_cache = (char **)malloc(sizeof(char *)*CACHESIZE);
132 lru = (int *)malloc(sizeof(int)*CACHESIZE);
133 for(k=0; k<CACHESIZE; k++)
135 data_cache[k] = NULL;
136 fname_cache[k] = NULL;
137 lru[k] = 0;
141 least = 0;
142 least_idx = 0;
143 for(k=0; k<CACHESIZE; k++)
145 lru[k]++;
146 if (lru[k] > least)
148 least = lru[k];
149 least_idx = k;
153 if (data_cache[least_idx] == NULL)
155 data_cache[least_idx] = localptr;
156 fname_cache[least_idx] = fname;
157 lru[least_idx] = 0;
159 else
161 free_newtile(data_cache[least_idx]);
162 data_cache[least_idx] = localptr;
163 free(fname_cache[least_idx]);
164 fname_cache[least_idx] = fname;
165 lru[least_idx] = 0;
168 retval = localptr;
169 return retval;
172 void build_supertile(int i, int j)
174 int ii, jj, kk;
175 int doflip;
176 int *** newtile;
178 if (i < 0)
179 i = i + IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON);
180 if (i >= IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON))
181 i = i - IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON);
182 if (j < 0)
183 j = j + IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT);
184 if (j >= IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT))
185 j = j - IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT);
187 if (supertile == NULL)
189 supertile = (int ***)malloc(sizeof(int **)*3*IN_TILE_PTS_X);
190 for(ii=0; ii<3*IN_TILE_PTS_X; ii++)
192 supertile[ii] = (int **)malloc(sizeof(int *)*3*IN_TILE_PTS_Y);
193 for(jj=0; jj<3*IN_TILE_PTS_Y; jj++)
195 supertile[ii][jj] = (int *)malloc(sizeof(int)*ZDIM);
200 supertile_min_x = (i / IN_TILE_PTS_X)*IN_TILE_PTS_X;
201 supertile_min_y = (j / IN_TILE_PTS_Y)*IN_TILE_PTS_Y;
203 /* Get tile containing (i,j) from cache*/
204 /* Get surrounding tiles from cache */
206 /* Lower-left */
207 ii = i - IN_TILE_PTS_X;
208 if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
209 doflip = 0;
210 jj = j - IN_TILE_PTS_Y;
211 if (jj < 0)
213 doflip = 1;
214 jj = j;
215 ii -= (int)(180./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
216 if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
218 newtile = get_tile_from_cache(ii, jj);
219 if (doflip)
221 for(ii=0; ii<IN_TILE_PTS_X; ii++)
223 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
225 for(kk=0; kk<ZDIM; kk++)
226 supertile[ii][jj][kk] = newtile[ii][IN_TILE_PTS_Y-1-jj][kk];
230 else
232 for(ii=0; ii<IN_TILE_PTS_X; ii++)
234 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
236 for(kk=0; kk<ZDIM; kk++)
237 supertile[ii][jj][kk] = newtile[ii][jj][kk];
242 /* Left */
243 ii = i - IN_TILE_PTS_X;
244 if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
245 jj = j;
246 newtile = get_tile_from_cache(ii, jj);
247 for(ii=0; ii<IN_TILE_PTS_X; ii++)
249 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
251 for(kk=0; kk<ZDIM; kk++)
252 supertile[ii][jj+IN_TILE_PTS_Y][kk] = newtile[ii][jj][kk];
256 /* Upper-left */
257 ii = i - IN_TILE_PTS_X;
258 if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
259 doflip = 0;
260 jj = j + IN_TILE_PTS_Y;
261 if (jj >= (int)(180./IN_TILE_DEGREES_LAT)*IN_TILE_PTS_Y)
263 doflip = 1;
264 jj = j;
265 ii -= (int)(180./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
266 if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
268 newtile = get_tile_from_cache(ii, jj);
269 if (doflip)
271 for(ii=0; ii<IN_TILE_PTS_X; ii++)
273 for(jj=0; jj<IN_TILE_PTS_X; jj++)
275 for(kk=0; kk<ZDIM; kk++)
276 supertile[ii][jj+2*IN_TILE_PTS_Y][kk] = newtile[ii][IN_TILE_PTS_Y-1-jj][kk];
280 else
282 for(ii=0; ii<IN_TILE_PTS_X; ii++)
284 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
286 for(kk=0; kk<ZDIM; kk++)
287 supertile[ii][jj+2*IN_TILE_PTS_Y][kk] = newtile[ii][jj][kk];
292 /* Below */
293 ii = i;
294 doflip = 0;
295 jj = j - IN_TILE_PTS_Y;
296 if (jj < 0)
298 doflip = 1;
299 jj = j;
300 ii -= (int)(180./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
301 if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
303 newtile = get_tile_from_cache(ii, jj);
304 if (doflip)
306 for(ii=0; ii<IN_TILE_PTS_X; ii++)
308 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
310 for(kk=0; kk<ZDIM; kk++)
311 supertile[ii+IN_TILE_PTS_X][jj][kk] = newtile[ii][IN_TILE_PTS_Y-1-jj][kk];
315 else
317 for(ii=0; ii<IN_TILE_PTS_X; ii++)
319 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
321 for(kk=0; kk<ZDIM; kk++)
322 supertile[ii+IN_TILE_PTS_X][jj][kk] = newtile[ii][jj][kk];
327 /* Center */
328 newtile = get_tile_from_cache(i, j);
329 for(ii=0; ii<IN_TILE_PTS_X; ii++)
331 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
333 for(kk=0; kk<ZDIM; kk++)
334 supertile[ii+IN_TILE_PTS_X][jj+IN_TILE_PTS_Y][kk] = newtile[ii][jj][kk];
338 /* Above */
339 ii = i;
340 doflip = 0;
341 jj = j + IN_TILE_PTS_Y;
342 if (jj >= (int)(180./IN_TILE_DEGREES_LAT)*IN_TILE_PTS_Y)
344 doflip = 1;
345 jj = j;
346 ii -= (int)(180./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
347 if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
349 newtile = get_tile_from_cache(ii, jj);
350 if (doflip)
352 for(ii=0; ii<IN_TILE_PTS_X; ii++)
354 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
356 for(kk=0; kk<ZDIM; kk++)
357 supertile[ii+IN_TILE_PTS_X][jj+2*IN_TILE_PTS_Y][kk] = newtile[ii][IN_TILE_PTS_Y-1-jj][kk];
361 else
363 for(ii=0; ii<IN_TILE_PTS_X; ii++)
365 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
367 for(kk=0; kk<ZDIM; kk++)
368 supertile[ii+IN_TILE_PTS_X][jj+2*IN_TILE_PTS_Y][kk] = newtile[ii][jj][kk];
373 /* Lower-right */
374 ii = i + IN_TILE_PTS_X;
375 if (ii >= (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X) ii -= (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
376 doflip = 0;
377 jj = j - IN_TILE_PTS_Y;
378 if (jj < 0)
380 doflip = 1;
381 jj = j;
382 ii -= (int)(180./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
383 if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
385 newtile = get_tile_from_cache(ii, jj);
386 if (doflip)
388 for(ii=0; ii<IN_TILE_PTS_X; ii++)
390 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
392 for(kk=0; kk<ZDIM; kk++)
393 supertile[ii+2*IN_TILE_PTS_X][jj][kk] = newtile[ii][IN_TILE_PTS_Y-1-jj][kk];
397 else
399 for(ii=0; ii<IN_TILE_PTS_X; ii++)
401 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
403 for(kk=0; kk<ZDIM; kk++)
404 supertile[ii+2*IN_TILE_PTS_X][jj][kk] = newtile[ii][jj][kk];
409 /* Right */
410 ii = i + IN_TILE_PTS_X;
411 if (ii >= (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X) ii -= (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
412 jj = j;
413 newtile = get_tile_from_cache(ii, jj);
414 for(ii=0; ii<IN_TILE_PTS_X; ii++)
416 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
418 for(kk=0; kk<ZDIM; kk++)
419 supertile[ii+2*IN_TILE_PTS_X][jj+IN_TILE_PTS_Y][kk] = newtile[ii][jj][kk];
423 /* Upper-right */
424 ii = i + IN_TILE_PTS_X;
425 if (ii >= (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X) ii -= (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
426 doflip = 0;
427 jj = j + IN_TILE_PTS_Y;
428 if (jj >= (int)(180./IN_TILE_DEGREES_LAT)*IN_TILE_PTS_Y)
430 doflip = 1;
431 jj = j;
432 ii -= (int)(180./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
433 if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
435 newtile = get_tile_from_cache(ii, jj);
436 if (doflip)
438 for(ii=0; ii<IN_TILE_PTS_X; ii++)
440 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
442 for(kk=0; kk<ZDIM; kk++)
443 supertile[ii+2*IN_TILE_PTS_X][jj+2*IN_TILE_PTS_Y][kk] = newtile[ii][IN_TILE_PTS_Y-1-jj][kk];
447 else
449 for(ii=0; ii<IN_TILE_PTS_X; ii++)
451 for(jj=0; jj<IN_TILE_PTS_Y; jj++)
453 for(kk=0; kk<ZDIM; kk++)
454 supertile[ii+2*IN_TILE_PTS_X][jj+2*IN_TILE_PTS_Y][kk] = newtile[ii][jj][kk];
460 int get_value(int i, int j, int k, int irad)
462 int i_src, j_src;
463 int ii, jj;
464 int n, sum;
465 float r_rel;
467 if (i < 0)
468 i = i + IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON);
469 if (i >= IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON))
470 i = i - IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON);
471 if (j < 0)
472 j = j + IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT);
473 if (j >= IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT))
474 j = j - IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT);
476 i_src = i % IN_TILE_PTS_X + IN_TILE_PTS_X;
477 j_src = j % IN_TILE_PTS_Y + IN_TILE_PTS_Y;
479 /* Interpolate values from supertile */
480 sum = 0;
481 n = 0;
482 for(ii=i_src; ii<=i_src+irad; ii++)
484 for(jj=j_src; jj<=j_src+irad; jj++)
486 sum += supertile[ii][jj][k];
487 n++;
491 if (n > 0) return sum/n;
492 else return 0;
495 int is_in_supertile(int i, int j)
497 if (i < 0)
498 i = i + IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON);
499 if (i >= IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON))
500 i = i - IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON);
501 if (j < 0)
502 j = j + IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT);
503 if (j >= IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT))
504 j = j - IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT);
506 /* Check whether (i,j) is in the interior of supertile */
507 if ((i >= supertile_min_x) && (i < supertile_min_x+IN_TILE_PTS_X) &&
508 (j >= supertile_min_y) && (j < supertile_min_y+IN_TILE_PTS_Y))
509 return 1;
510 else
511 return 0;
514 int main(int argc, char ** argv)
516 int tile_x, tile_y, input_x, input_y, temp_x, temp_y;
517 int i, j, k, z, ii, jj;
518 int i_src, j_src;
519 int *** intdata;
520 int out_fd;
521 int ir_rel;
522 float r_rel;
523 unsigned char * outdata;
524 char out_filename[256];
526 r_rel = (float)IN_TILE_PTS_X/(float)OUT_TILE_PTS_X*OUT_TILE_DEGREES_LON/IN_TILE_DEGREES_LON;
527 ir_rel = (int)rint(r_rel);
529 /* Allocate memory to hold a single output tile */
530 intdata = (int ***)malloc(sizeof(int **)*(OUT_TILE_PTS_X+2*HALO_WIDTH));
531 for(i=0; i<(OUT_TILE_PTS_X+2*HALO_WIDTH); i++)
533 intdata[i] = (int **)malloc(sizeof(int *)*(OUT_TILE_PTS_Y+2*HALO_WIDTH));
534 for(j=0; j<(OUT_TILE_PTS_Y+2*HALO_WIDTH); j++)
536 intdata[i][j] = (int *)malloc(sizeof(int)*ZDIM);
540 /* Allocate output buffer */
541 outdata = (unsigned char *)malloc((OUT_TILE_PTS_X+2*HALO_WIDTH)*(OUT_TILE_PTS_Y+2*HALO_WIDTH)*ZDIM*N_BYTES);
543 for(tile_x=0; tile_x<OUT_TILE_PTS_X*(int)(360.0/OUT_TILE_DEGREES_LON); tile_x+=OUT_TILE_PTS_X)
545 for(tile_y=0; tile_y<OUT_TILE_PTS_Y*(int)(180.0/OUT_TILE_DEGREES_LAT); tile_y+=OUT_TILE_PTS_Y)
547 /* Build name of output file for current tile_x and tile_y */
548 sprintf(out_filename,"retiled/%5.5i-%5.5i.%5.5i-%5.5i",tile_x+1,tile_x+OUT_TILE_PTS_X,tile_y+1,tile_y+OUT_TILE_PTS_Y);
550 /* Initialize the output data for current tile */
551 for(i=0; i<(OUT_TILE_PTS_X+2*HALO_WIDTH); i++)
553 for(j=0; j<(OUT_TILE_PTS_Y+2*HALO_WIDTH); j++)
555 intdata[i][j][0] = MSG_FLAG;
559 /* Attempt to open output file */
560 if ((out_fd = open(out_filename, O_WRONLY|O_CREAT|O_TRUNC, S_IRUSR|S_IWUSR|S_IRGRP|S_IROTH)) == -1)
562 fprintf(stderr, "Error: Could not create or open file %s\n", out_filename);
563 return 1;
566 /* Fill tile with data */
567 for(i=-1*HALO_WIDTH; i<OUT_TILE_PTS_X+HALO_WIDTH; i++)
569 for(j=-1*HALO_WIDTH; j<OUT_TILE_PTS_Y+HALO_WIDTH; j++)
571 if (intdata[i+HALO_WIDTH][j+HALO_WIDTH][0] == MSG_FLAG)
573 i_src = ir_rel*(tile_x+i);
574 j_src = ir_rel*(tile_y+j);
576 build_supertile(i_src,j_src);
578 for(ii=-1*HALO_WIDTH; ii<OUT_TILE_PTS_X+HALO_WIDTH; ii++)
580 for(jj=-1*HALO_WIDTH; jj<OUT_TILE_PTS_Y+HALO_WIDTH; jj++)
582 i_src = ir_rel*(tile_x+ii);
583 j_src = ir_rel*(tile_y+jj);
584 if (is_in_supertile(i_src,j_src))
586 for(k=0; k<ZDIM; k++)
587 intdata[ii+HALO_WIDTH][jj+HALO_WIDTH][k] = get_value(i_src,j_src,k,ir_rel-1);
596 /* Write out the data */
597 for(i=0; i<(OUT_TILE_PTS_X+2*HALO_WIDTH); i++)
599 for(j=0; j<(OUT_TILE_PTS_Y+2*HALO_WIDTH); j++)
601 for(z=0; z<ZDIM; z++)
603 for(k=0; k<N_BYTES; k++)
605 outdata[z*N_BYTES*(OUT_TILE_PTS_Y+2*HALO_WIDTH)*(OUT_TILE_PTS_X+2*HALO_WIDTH)+j*N_BYTES*(OUT_TILE_PTS_X+2*HALO_WIDTH)+i*N_BYTES+k] =
606 (intdata[i][j][z] & (0xff << 8*(N_BYTES-k-1))) >> (8*(N_BYTES-k-1));
611 write(out_fd,(void *)outdata,(OUT_TILE_PTS_X+2*HALO_WIDTH)*(OUT_TILE_PTS_Y+2*HALO_WIDTH)*ZDIM*N_BYTES);
612 close(out_fd);
613 printf("Wrote file %s\n",out_filename);
617 /* Deallocate memory */
618 for(i=0; i<(OUT_TILE_PTS_X+2*HALO_WIDTH); i++)
620 for(j=0; j<(OUT_TILE_PTS_Y+2*HALO_WIDTH); j++)
622 free(intdata[i][j]);
624 free(intdata[i]);
626 free(intdata);
629 return 0;