Skip to content
GitLab
Commit 70f19ca4 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

remapbil, remapbic: improvement and speedup for regional lonlat grids

parent b0caeaca
Hide whitespace changes
Inline Side-by-side
ChangeLog
View file @ 70f19ca4
......@@ -2,6 +2,7 @@
2006-09-?? Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* using CDI library version 1.0.2
* remapbil, remapbic: speedup for regional lonlat grids
* Seltime: print warning message if parameter not found
* pstreamDefVlist: unpack netCDF data if datatype = FLT64
* Version 1.0.2 released
......
config/default
View file @ 70f19ca4
......@@ -97,19 +97,19 @@ case "${HOSTNAME}" in
ecgate)
./configure --prefix=$HOME/local \
--with-netcdf=/home/ms/spdekplb/hmk/local/AIX \
CC=xlc_r CFLAGS="-g -O -q64 -qMAXMEM=-1 -qarch=auto -DHAVE_MMAP -D_LARGE_FILES=1"
CC=xlc_r CFLAGS="-O -q64 -qMAXMEM=-1 -qarch=auto -DHAVE_MMAP -D_LARGE_FILES=1"
;;
# powerpc-ibm-aix5.1.0.0
hpcd*)
# powerpc-ibm-aix5.3.0.0
hpc*)
./configure --prefix=$HOME/local \
--with-netcdf=/usr/local/lib/netcdf-3.5.0/LP64 \
CC=xlc_r CFLAGS="-g -O -q64 -qMAXMEM=-1 -qarch=auto -DHAVE_MMAP -D_LARGE_FILES=1"
CC=xlc_r CFLAGS="-O3 -q64 -qstrict -qarch=auto -qtune=auto -DHAVE_MMAP -D_LARGE_FILES=1"
;;
# powerpc-ibm-aix5.2.0.0
p020et01)
./configure --prefix=$HOME/local \
--with-netcdf=/uhome/mpikornb \
CC=xlc_r CFLAGS="-g -O -q64 -qMAXMEM=-1 -qarch=auto -DHAVE_MMAP -D_LARGE_FILES=1"
CC=xlc_r CFLAGS="-O -q64 -qMAXMEM=-1 -qarch=auto -DHAVE_MMAP -D_LARGE_FILES=1"
;;
# hppa2.0w-hp-hpux11.11
moon)
......
src/Remap.c
View file @ 70f19ca4
......@@ -66,6 +66,7 @@ void *Remap(void *argument)
int max_remaps = 0;
int remap_test = 0;
int lgridboxinfo = TRUE;
int grid1sizemax;
char varname[128];
double missval;
double *array1 = NULL, *array2 = NULL;
......@@ -187,6 +188,8 @@ void *Remap(void *argument)
vlistChangeGridIndex(vlistID2, index, gridID2);
}
gridID1 = vlistGrid(vlistID1, 0);
if ( max_remaps == 0 )
{
nzaxis = vlistNzaxis(vlistID1);
......@@ -219,7 +222,13 @@ void *Remap(void *argument)
remaps[0].gridID = gridID1;
remaps[0].gridsize = gridInqSize(gridID1);
remaps[0].nmiss = 0;
if ( gridInqType(gridID1) == GRID_LONLAT && gridIsRotated(gridID1) &&
map_type != MAP_TYPE_CONSERV )
remaps[0].gridsize += 4*(gridInqXsize(gridID1)+2) + 4*(gridInqYsize(gridID1)+2);
if ( gridInqType(gridID1) == GRID_GME ) gridsize = remaps[0].grid.grid1_nvgp;
if ( gridsize != remaps[0].gridsize )
cdoAbort("Size of data grid and weights from %s differ!", remap_file);
......@@ -291,17 +300,17 @@ void *Remap(void *argument)
cdoAbort("unknown mapping method");
}
gridsize = vlistGridsizeMax(vlistID1);
grid1sizemax = vlistGridsizeMax(vlistID1);
if ( map_type == MAP_TYPE_BICUBIC )
{
grad1_lat = (double *) malloc(gridsize*sizeof(double));
grad1_lon = (double *) malloc(gridsize*sizeof(double));
grad1_latlon = (double *) malloc(gridsize*sizeof(double));
grad1_lat = (double *) malloc(grid1sizemax*sizeof(double));
grad1_lon = (double *) malloc(grid1sizemax*sizeof(double));
grad1_latlon = (double *) malloc(grid1sizemax*sizeof(double));
}
array1 = (double *) malloc(gridsize*sizeof(double));
imask = (int *) malloc(gridsize*sizeof(int));
array1 = (double *) malloc(grid1sizemax*sizeof(double));
imask = (int *) malloc(grid1sizemax*sizeof(int));
gridsize = gridInqSize(gridID2);
array2 = (double *) malloc(gridsize*sizeof(double));
......@@ -339,6 +348,56 @@ void *Remap(void *argument)
missval = vlistInqVarMissval(vlistID1, varID);
gridsize = gridInqSize(gridID1);
if ( gridInqType(gridID1) == GRID_LONLAT && gridIsRotated(gridID1) &&
map_type != MAP_TYPE_CONSERV )
{
int gridsize_new;
int nx, ny;
nx = gridInqXsize(gridID1);
ny = gridInqYsize(gridID1);
/* gridsize_new = gridsize + 2*(nx+1) + 2*(ny+1); */
gridsize_new = gridsize + 4*(nx+2) + 4*(ny+2);
/*
fprintf(stderr, "Data on rotated grid found!\n");
fprintf(stderr, "gridsize %d %d %d %d\n", gridsize, gridsize_new, nx, ny);
*/
if ( gridsize_new > grid1sizemax )
{
grid1sizemax = gridsize_new;
array1 = (double *) realloc(array1, grid1sizemax*sizeof(double));
imask = (int *) realloc(imask, grid1sizemax*sizeof(int));
if ( map_type == MAP_TYPE_BICUBIC )
{
grad1_lat = (double *) realloc(grad1_lat, grid1sizemax*sizeof(double));
grad1_lon = (double *) realloc(grad1_lon, grid1sizemax*sizeof(double));
grad1_latlon = (double *) realloc(grad1_latlon, grid1sizemax*sizeof(double));
}
}
for ( j = ny-1; j >= 0; j-- )
for ( i = nx-1; i >= 0; i-- )
array1[(j+2)*(nx+4)+i+2] = array1[j*nx+i];
/* array1[(j+1)*(nx+2)+i+1] = array1[j*nx+i];*/
/*
for ( j = 0; j < ny+2; j++ ) array1[j*(nx+2)+0] = missval;
for ( j = 0; j < ny+2; j++ ) array1[j*(nx+2)+nx+1] = missval;
for ( i = 0; i < nx+2; i++ ) array1[ 0*(nx+2)+i] = missval;
for ( i = 0; i < nx+2; i++ ) array1[(ny+1)*(nx+2)+i] = missval;
*/
for ( j = 0; j < ny+4; j++ ) array1[j*(nx+4)+0] = missval;
for ( j = 0; j < ny+4; j++ ) array1[j*(nx+4)+1] = missval;
for ( j = 0; j < ny+4; j++ ) array1[j*(nx+4)+nx+2] = missval;
for ( j = 0; j < ny+4; j++ ) array1[j*(nx+4)+nx+3] = missval;
for ( i = 0; i < nx+4; i++ ) array1[ 0*(nx+4)+i] = missval;
for ( i = 0; i < nx+4; i++ ) array1[ 1*(nx+4)+i] = missval;
for ( i = 0; i < nx+4; i++ ) array1[(ny+2)*(nx+4)+i] = missval;
for ( i = 0; i < nx+4; i++ ) array1[(ny+3)*(nx+4)+i] = missval;
gridsize = gridsize_new;
nmiss1 += 4*(nx+2) + 4*(ny+2);
}
for ( i = 0; i < gridsize; i++ )
if ( DBL_IS_EQUAL(array1[i], missval) )
imask[i] = FALSE;
......
src/remap.h
View file @ 70f19ca4
#define NORM_OPT_NONE 1
#define NORM_OPT_DESTAREA 2
#define NORM_OPT_FRACAREA 3
......@@ -21,6 +22,7 @@ typedef struct {
int pinit; /* TRUE if the pointers are initialized */
int gridID1;
int gridID2;
int no_fall_back;
int grid1_size, grid2_size; /* total points on each grid */
int grid1_rank, grid2_rank; /* rank of each grid */
int grid1_corners, grid2_corners; /* number of corners for each grid cell */
......
src/remaplib.c
View file @ 70f19ca4
......@@ -362,6 +362,8 @@ void remapGridInit(int map_type, int gridID1, int gridID2, REMAPGRID *rg)
double tmp_lats[4], tmp_lons[4]; /* temps for computing bounding boxes */
rg->no_fall_back = FALSE;
if ( map_type == MAP_TYPE_CONSERV )
{
rg->luse_grid1_corners = TRUE;
......@@ -391,6 +393,85 @@ void remapGridInit(int map_type, int gridID1, int gridID2, REMAPGRID *rg)
}
else
{
if ( gridInqType(gridID1) == GRID_LONLAT && gridIsRotated(gridID1) &&
map_type != MAP_TYPE_CONSERV )
{
/*
int gridIDnew;
int nx, ny, nxp2, nyp2;
double *xvals, *yvals;
nx = gridInqXsize(gridID1);
ny = gridInqYsize(gridID1);
nxp2 = nx+2;
nyp2 = ny+2;
xvals = (double *) malloc(nxp2*sizeof(double));
yvals = (double *) malloc(nyp2*sizeof(double));
gridInqXvals(gridID1, xvals+1);
gridInqYvals(gridID1, yvals+1);
gridIDnew = gridCreate(GRID_LONLAT, nxp2*nyp2);
gridDefXsize(gridIDnew, nxp2);
gridDefYsize(gridIDnew, nyp2);
xvals[0] = xvals[1] - gridInqXinc(gridID1);
xvals[nxp2-1] = xvals[nx] + gridInqXinc(gridID1);
yvals[0] = yvals[1] - gridInqYinc(gridID1);
yvals[nyp2-1] = yvals[ny] + gridInqYinc(gridID1);
gridDefXvals(gridIDnew, xvals);
gridDefYvals(gridIDnew, yvals);
free(xvals);
free(yvals);
gridDefXpole(gridIDnew, gridInqXpole(gridID1));
gridDefYpole(gridIDnew, gridInqYpole(gridID1));
*/
int gridIDnew;
int nx, ny, nxp4, nyp4;
double *xvals, *yvals;
nx = gridInqXsize(gridID1);
ny = gridInqYsize(gridID1);
nxp4 = nx+4;
nyp4 = ny+4;
xvals = (double *) malloc(nxp4*sizeof(double));
yvals = (double *) malloc(nyp4*sizeof(double));
gridInqXvals(gridID1, xvals+2);
gridInqYvals(gridID1, yvals+2);
gridIDnew = gridCreate(GRID_LONLAT, nxp4*nyp4);
gridDefXsize(gridIDnew, nxp4);
gridDefYsize(gridIDnew, nyp4);
xvals[0] = xvals[2] - 2*gridInqXinc(gridID1);
xvals[1] = xvals[2] - gridInqXinc(gridID1);
xvals[nxp4-2] = xvals[nx+1] + gridInqXinc(gridID1);
xvals[nxp4-1] = xvals[nx+1] + 2*gridInqXinc(gridID1);
yvals[0] = yvals[2] - 2*gridInqYinc(gridID1);
yvals[1] = yvals[2] - gridInqYinc(gridID1);
yvals[nyp4-2] = yvals[ny+1] + gridInqYinc(gridID1);
yvals[nyp4-1] = yvals[ny+1] + 2*gridInqYinc(gridID1);
gridDefXvals(gridIDnew, xvals);
gridDefYvals(gridIDnew, yvals);
free(xvals);
free(yvals);
gridDefXpole(gridIDnew, gridInqXpole(gridID1));
gridDefYpole(gridIDnew, gridInqYpole(gridID1));
rg->gridID1 = gridIDnew;
rg->no_fall_back = TRUE;
}
gridID1 = gridToCurvilinear(rg->gridID1);
lgrid1_gen_bounds = TRUE;
}
......@@ -1466,6 +1547,8 @@ void grid_search(REMAPGRID *rg, int *src_add, double *src_lats, double *src_lons
bilinear weights.
*/
if ( rg->no_fall_back ) return;
/*
printf("Could not find location for %g %g\n", plat*RAD2DEG, plon*RAD2DEG);
printf("Using nearest-neighbor average for this point\n");
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Imprint | Privacy Policy