Commit 52e210e6 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

gridToUnstructured: optional computation of grid bounds

parent abd4e77b
......@@ -96,7 +96,7 @@ void *Arithlat(void *argument)
gridtype == GRID_GAUSSIAN ||
gridtype == GRID_LCC )
{
gridID = gridToCurvilinear(gridID, 1);
gridID = gridToCurvilinear(gridID, 0);
}
else if ( gridtype == GRID_CURVILINEAR ||
gridtype == GRID_UNSTRUCTURED )
......@@ -105,7 +105,7 @@ void *Arithlat(void *argument)
}
else if ( gridtype == GRID_GME )
{
gridID = gridToUnstructured(gridID);
gridID = gridToUnstructured(gridID, 0);
}
else
{
......
......@@ -155,7 +155,7 @@ void *Output(void *argument)
if ( operatorID == OUTPUTFLD || operatorID == OUTPUTXYZ || operatorID == OUTPUTTAB )
{
if ( gridInqType(gridID) == GRID_GME ) gridID = gridToUnstructured(gridID);
if ( gridInqType(gridID) == GRID_GME ) gridID = gridToUnstructured(gridID, 0);
if ( gridInqType(gridID) != GRID_UNSTRUCTURED && gridInqType(gridID) != GRID_CURVILINEAR )
gridID = gridToCurvilinear(gridID, 0);
......
......@@ -656,7 +656,7 @@ void *Outputgmt(void *argument)
zaxisID = vlistInqVarZaxis(vlistID, varID);
missval = vlistInqVarMissval(vlistID, varID);
if ( gridInqType(gridID) == GRID_GME ) gridID = gridToUnstructured(gridID);
if ( gridInqType(gridID) == GRID_GME ) gridID = gridToUnstructured(gridID, 1);
if ( gridInqType(gridID) != GRID_UNSTRUCTURED && gridInqType(gridID) != GRID_CURVILINEAR )
{
......
......@@ -555,7 +555,7 @@ void *Remap(void *argument)
remaps[0].grid.grid2_nvgp = gridInqSize(gridID2);
remaps[0].grid.grid2_vgpm = (int *) realloc(remaps[0].grid.grid2_vgpm,
gridInqSize(gridID2)*sizeof(int));
gridID2_gme = gridToUnstructured(gridID2);
gridID2_gme = gridToUnstructured(gridID2, 1);
gridInqMaskGME(gridID2_gme, remaps[0].grid.grid2_vgpm);
for ( i = 0; i < gridsize2; ++i )
if ( remaps[0].grid.grid2_vgpm[i] ) isize++;
......
......@@ -244,7 +244,7 @@ void *Setgrid(void *argument)
else if ( gridtype == GRID_UNSTRUCTURED )
{
if ( gridInqType(gridID1) == GRID_GME ) ligme = 1;
gridID2 = gridToUnstructured(gridID1);
gridID2 = gridToUnstructured(gridID1, 1);
if ( ligme )
{
......
......@@ -47,7 +47,7 @@ void *Writegrid(void *argument)
gridtype = gridInqType(gridID);
gridsize = gridInqSize(gridID);
if ( gridtype == GRID_GME ) gridID = gridToUnstructured(gridID);
if ( gridtype == GRID_GME ) gridID = gridToUnstructured(gridID, 1);
if ( gridtype != GRID_CURVILINEAR && gridtype != GRID_UNSTRUCTURED )
gridID = gridToCurvilinear(gridID, 1);
......
......@@ -1050,7 +1050,7 @@ void eca4(const ECA_REQUEST_4 *request)
if ( gridtype != GRID_UNSTRUCTURED && gridtype != GRID_CURVILINEAR )
{
if ( gridtype == GRID_GME )
gridID = gridToUnstructured(gridID);
gridID = gridToUnstructured(gridID, 1);
else
gridID = gridToCurvilinear(gridID, 1);
}
......
......@@ -853,7 +853,7 @@ int gridToCurvilinear(int gridID1, int lbounds)
}
int gridToUnstructured(int gridID1)
int gridToUnstructured(int gridID1, int lbounds)
{
int gridID2;
int gridtype, gridsize;
......@@ -872,8 +872,6 @@ int gridToUnstructured(int gridID1)
int nx, ny;
double *xvals, *yvals;
double *xvals2D, *yvals2D;
double *xbounds = NULL, *ybounds = NULL;
double *xbounds2D, *ybounds2D;
gridDefXname(gridID2, "lon");
gridDefYname(gridID2, "lat");
......@@ -930,54 +928,60 @@ int gridToUnstructured(int gridID1)
free(xvals2D);
free(yvals2D);
if ( gridInqXbounds(gridID1, NULL) )
if ( lbounds )
{
xbounds = (double *) malloc(2*nx*sizeof(double));
gridInqXbounds(gridID1, xbounds);
}
else if ( nx > 1 )
{
xbounds = (double *) malloc(2*nx*sizeof(double));
gridGenXbounds(nx, xvals, xbounds);
}
double *xbounds = NULL, *ybounds = NULL;
double *xbounds2D, *ybounds2D;
if ( gridInqYbounds(gridID1, NULL) )
{
ybounds = (double *) malloc(2*ny*sizeof(double));
gridInqYbounds(gridID1, ybounds);
}
else if ( ny > 1 )
{
ybounds = (double *) malloc(2*ny*sizeof(double));
gridGenYbounds(ny, yvals, ybounds);
}
free(xvals);
free(yvals);
if ( xbounds && ybounds )
{
xbounds2D = (double *) malloc(4*gridsize*sizeof(double));
ybounds2D = (double *) malloc(4*gridsize*sizeof(double));
if ( gridInqXbounds(gridID1, NULL) )
{
xbounds = (double *) malloc(2*nx*sizeof(double));
gridInqXbounds(gridID1, xbounds);
}
else if ( nx > 1 )
{
xbounds = (double *) malloc(2*nx*sizeof(double));
gridGenXbounds(nx, xvals, xbounds);
}
if ( gridIsRotated(gridID1) )
if ( gridInqYbounds(gridID1, NULL) )
{
gridGenRotBounds(gridID1, nx, ny, xbounds, ybounds, xbounds2D, ybounds2D);
ybounds = (double *) malloc(2*ny*sizeof(double));
gridInqYbounds(gridID1, ybounds);
}
else
else if ( ny > 1 )
{
gridGenXbounds2D(nx, ny, xbounds, xbounds2D);
gridGenYbounds2D(nx, ny, ybounds, ybounds2D);
ybounds = (double *) malloc(2*ny*sizeof(double));
gridGenYbounds(ny, yvals, ybounds);
}
gridDefXbounds(gridID2, xbounds2D);
gridDefYbounds(gridID2, ybounds2D);
if ( xbounds && ybounds )
{
xbounds2D = (double *) malloc(4*gridsize*sizeof(double));
ybounds2D = (double *) malloc(4*gridsize*sizeof(double));
free(xbounds);
free(ybounds);
free(xbounds2D);
free(ybounds2D);
}
if ( gridIsRotated(gridID1) )
{
gridGenRotBounds(gridID1, nx, ny, xbounds, ybounds, xbounds2D, ybounds2D);
}
else
{
gridGenXbounds2D(nx, ny, xbounds, xbounds2D);
gridGenYbounds2D(nx, ny, ybounds, ybounds2D);
}
gridDefXbounds(gridID2, xbounds2D);
gridDefYbounds(gridID2, ybounds2D);
free(xbounds);
free(ybounds);
free(xbounds2D);
free(ybounds2D);
}
}
free(xvals);
free(yvals);
gridCopyMask(gridID1, gridID2, gridsize);
......@@ -999,7 +1003,7 @@ int gridToUnstructured(int gridID1)
int nv = 6;
int *imask;
double *xvals, *yvals;
double *xbounds, *ybounds;
double *xbounds = NULL, *ybounds = NULL;
nd = gridInqGMEnd(gridID1);
ni = gridInqGMEni(gridID1);
......@@ -1009,21 +1013,25 @@ int gridToUnstructured(int gridID1)
imask = (int *) malloc(gridsize*sizeof(int));
xvals = (double *) malloc(gridsize*sizeof(double));
yvals = (double *) malloc(gridsize*sizeof(double));
xbounds = (double *) malloc(nv*gridsize*sizeof(double));
ybounds = (double *) malloc(nv*gridsize*sizeof(double));
if ( lbounds )
{
xbounds = (double *) malloc(nv*gridsize*sizeof(double));
ybounds = (double *) malloc(nv*gridsize*sizeof(double));
}
gme_grid(gridsize, xvals, yvals, xbounds, ybounds, imask, ni, nd, ni2, ni3);
gme_grid(lbounds, gridsize, xvals, yvals, xbounds, ybounds, imask, ni, nd, ni2, ni3);
for ( i = 0; i < gridsize; i++ )
{
xvals[i] *= RAD2DEG;
yvals[i] *= RAD2DEG;
for ( j = 0; j < nv; j++ )
{
xbounds[i*nv + j] *= RAD2DEG;
ybounds[i*nv + j] *= RAD2DEG;
}
if ( lbounds )
for ( j = 0; j < nv; j++ )
{
xbounds[i*nv + j] *= RAD2DEG;
ybounds[i*nv + j] *= RAD2DEG;
}
/* printf("%d %g %g\n", i, xvals[i], yvals[i]); */
}
......@@ -1037,8 +1045,11 @@ int gridToUnstructured(int gridID1)
gridDefNvertex(gridID2, nv);
gridDefXbounds(gridID2, xbounds);
gridDefYbounds(gridID2, ybounds);
if ( lbounds )
{
gridDefXbounds(gridID2, xbounds);
gridDefYbounds(gridID2, ybounds);
}
gridDefXunits(gridID2, "degrees_east");
gridDefYunits(gridID2, "degrees_north");
......@@ -1046,8 +1057,8 @@ int gridToUnstructured(int gridID1)
free (imask);
free (xvals);
free (yvals);
free (xbounds);
free (ybounds);
if ( xbounds ) free (xbounds);
if ( ybounds ) free (ybounds);
gridCopyMask(gridID1, gridID2, gridsize);
......@@ -1390,7 +1401,7 @@ int gridGenArea(int gridID, double *area)
if ( gridtype == GRID_GME )
{
lgriddestroy = TRUE;
gridID = gridToUnstructured(gridID);
gridID = gridToUnstructured(gridID, 1);
grid_mask = (int *) malloc(gridsize*sizeof(int));
gridInqMaskGME(gridID, grid_mask);
}
......@@ -1517,7 +1528,7 @@ int gridGenWeights(int gridID, double *grid_area, double *grid_wgts)
if ( gridtype == GRID_GME )
{
gridID = gridToUnstructured(gridID);
gridID = gridToUnstructured(gridID, 1);
grid_mask = (int *) malloc(gridsize*sizeof(int));
gridInqMaskGME(gridID, grid_mask);
}
......
......@@ -9,7 +9,7 @@ int referenceToGrid(int gridID);
void gridToDegree(const char *units, const char *string, int gridsize, double *array);
int gridToZonal(int gridID);
int gridToMeridional(int gridID);
int gridToUnstructured(int gridID);
int gridToUnstructured(int gridID, int lbounds);
int gridToCurvilinear(int gridID, int lbounds);
int gridCurvilinearToRegular(int gridID);
int gridToRegular(int gridID);
......@@ -30,7 +30,7 @@ void correct_sinxvals(int xsize, int ysize, double *xvals);
struct cart gc2cc(struct geo *position);
void factorni(int kni, int *kni2, int *kni3);
void gme_grid_restore(double *p, int ni, int nd);
void gme_grid(int gridsize, double *rlon, double *rlat,
void gme_grid(int lbounds, int gridsize, double *rlon, double *rlat,
double *blon, double *blat, int *imask,
int ni, int nd, int ni2, int ni3);
......
......@@ -1319,7 +1319,7 @@ void gme_grid_restore(double *p, int ni, int nd)
/*****************************************************************************/
void gme_grid(int gridsize, double *rlon, double *rlat,
void gme_grid(int lbounds, int gridsize, double *rlon, double *rlat,
double *blon, double *blat, int *imask,
int ni, int nd, int ni2, int ni3)
{
......@@ -1327,7 +1327,6 @@ void gme_grid(int gridsize, double *rlon, double *rlat,
int i, j;
double *xn;
double *rlonx, *rlatx;
struct polygon *poly;
/* check gridsize */
if ( (ni+1)*(ni+1)*nd != gridsize )
......@@ -1344,7 +1343,6 @@ void gme_grid(int gridsize, double *rlon, double *rlat,
xn = (double *) malloc(gridsize*3*sizeof(double));
rlonx = (double *) malloc((ni+3)*(ni+3)*nd*sizeof(double));
rlatx = (double *) malloc((ni+3)*(ni+3)*nd*sizeof(double));
poly = (struct polygon *) malloc((ni+1)*(ni+1)*nd*sizeof(struct polygon));
im1s = 0;
im1e = ni;
......@@ -1358,26 +1356,33 @@ void gme_grid(int gridsize, double *rlon, double *rlat,
initmask(imask,ni,nd);
neighbours(rlonx,rlatx,im1s-1,im1e+1,im2s-1,im2e+1,nd,
poly,im1s,im1e,im2s,im2e,nd);
if ( lbounds )
{
struct polygon *poly;
poly = (struct polygon *) malloc((ni+1)*(ni+1)*nd*sizeof(struct polygon));
boundary(poly,im1s,im1e,im2s,im2e,nd);
neighbours(rlonx,rlatx,im1s-1,im1e+1,im2s-1,im2e+1,nd, poly,im1s,im1e,im2s,im2e,nd);
for ( i = 0; i < gridsize; i++ )
{
for ( j = 0; j < poly[i].type; j++ )
{
blon[i*6+j] = poly[i].boundary[j].lon;
blat[i*6+j] = poly[i].boundary[j].lat;
}
if ( poly[i].type == pentagon )
boundary(poly,im1s,im1e,im2s,im2e,nd);
for ( i = 0; i < gridsize; i++ )
{
blon[i*6+5] = blon[i*6+4];
blat[i*6+5] = blat[i*6+4];
for ( j = 0; j < poly[i].type; j++ )
{
blon[i*6+j] = poly[i].boundary[j].lon;
blat[i*6+j] = poly[i].boundary[j].lat;
}
if ( poly[i].type == pentagon )
{
blon[i*6+5] = blon[i*6+4];
blat[i*6+5] = blat[i*6+4];
}
}
free(poly);
}
free(poly);
free(rlatx);
free(rlonx);
free(xn);
......
......@@ -209,7 +209,7 @@ void intgrid(field_t *field1, field_t *field2)
int gridID1, gridID2;
int i, nmiss;
double *lon1, *lat1;
double **array1_2D;
double **array1_2D = NULL;
double **field;
double *array = NULL;
double *array1, *array2;
......@@ -235,6 +235,10 @@ void intgrid(field_t *field1, field_t *field2)
nlon2 = gridInqXsize(gridID2);
nlat2 = gridInqYsize(gridID2);
array1_2D = (double **) malloc(nlat1*sizeof(double *));
for ( ilat = 0; ilat < nlat1; ilat++ )
array1_2D[ilat] = array1 + ilat*nlon1;
if ( nlon2 == 1 && nlat2 == 1 )
{
double lon2, lat2;
......@@ -279,6 +283,8 @@ void intgrid(field_t *field1, field_t *field2)
int gridsize2;
double *lon2, *lat2;
if ( gridInqType(gridID2) == GRID_GME ) gridID2 = gridToUnstructured(gridID2, 0);
if ( gridInqType(gridID2) != GRID_UNSTRUCTURED && gridInqType(gridID2) != GRID_CURVILINEAR )
gridID2 = gridToCurvilinear(gridID2, 0);
......@@ -292,10 +298,6 @@ void intgrid(field_t *field1, field_t *field2)
gridInqXvals(gridID2, lon2);
gridInqYvals(gridID2, lat2);
array1_2D = (double **) malloc(nlat1*sizeof(double *));
for ( ilat = 0; ilat < nlat1; ilat++ )
array1_2D[ilat] = array1 + ilat*nlon1;
intlinarr2(missval,
nlon1, nlat1, array1_2D, lon1, lat1,
gridsize2, array2, lon2, lat2);
......
......@@ -1033,7 +1033,7 @@ void remapGridInit(int map_type, int lextrapolate, int gridID1, int gridID2, rem
{
if ( gridInqType(rg->gridID1) == GRID_GME )
{
gridID1_gme = gridToUnstructured(rg->gridID1);
gridID1_gme = gridToUnstructured(rg->gridID1, 1);
rg->grid1_nvgp = gridInqSize(gridID1_gme);
gridID1 = gridDuplicate(gridID1_gme);
gridCompress(gridID1);
......@@ -1051,7 +1051,7 @@ void remapGridInit(int map_type, int lextrapolate, int gridID1, int gridID2, rem
{
if ( gridInqType(rg->gridID2) == GRID_GME )
{
gridID2_gme = gridToUnstructured(rg->gridID2);
gridID2_gme = gridToUnstructured(rg->gridID2, 1);
rg->grid2_nvgp = gridInqSize(gridID2_gme);
gridID2 = gridDuplicate(gridID2_gme);
gridCompress(gridID2);
......@@ -7150,7 +7150,7 @@ void read_remap_scrip(const char *interp_file, int gridID1, int gridID2, int *ma
if ( gridInqType(gridID1) == GRID_GME )
{
rg->grid1_nvgp = gridInqSize(gridID1);
gridID1_gme_c = gridToUnstructured(gridID1);
gridID1_gme_c = gridToUnstructured(gridID1, 1);
}
remapGridRealloc(rv->map_type, rg);
......
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