Commit ed84c4c2 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

added support for curvilinear and unstructured predefined grid names

parent afb65f34
......@@ -163,9 +163,9 @@ int gridToMeridional(int gridID1)
}
void grid_gen_corners(long n, const double* restrict vals, double* restrict corners)
void grid_gen_corners(int n, const double* restrict vals, double* restrict corners)
{
long i;
int i;
if ( n == 1 )
{
......@@ -183,11 +183,9 @@ void grid_gen_corners(long n, const double* restrict vals, double* restrict corn
}
void grid_gen_bounds(long n, const double* restrict vals, double* restrict bounds)
void grid_gen_bounds(int n, const double *restrict vals, double *restrict bounds)
{
long i;
for ( i = 0; i < n-1; ++i )
for ( int i = 0; i < n-1; ++i )
{
bounds[2*i+1] = 0.5*(vals[i] + vals[i+1]);
bounds[2*(i+1)] = 0.5*(vals[i] + vals[i+1]);
......@@ -378,10 +376,10 @@ void gridGenRotBounds(int gridID, int nx, int ny,
}
}
static
void gridGenXbounds2D(long nx, long ny, const double* restrict xbounds, double* restrict xbounds2D)
void grid_gen_xbounds2D(int nx, int ny, const double *restrict xbounds, double *restrict xbounds2D)
{
long index;
int index;
double minlon, maxlon;
#if defined(_OPENMP)
......@@ -389,12 +387,12 @@ void gridGenXbounds2D(long nx, long ny, const double* restrict xbounds, double*
shared(nx, ny, xbounds, xbounds2D) \
private(minlon, maxlon, index)
#endif
for ( long i = 0; i < nx; ++i )
for ( int i = 0; i < nx; ++i )
{
minlon = xbounds[2*i ];
maxlon = xbounds[2*i+1];
for ( long j = 0; j < ny; ++j )
for ( int j = 0; j < ny; ++j )
{
index = j*4*nx + 4*i;
xbounds2D[index ] = minlon;
......@@ -405,10 +403,10 @@ void gridGenXbounds2D(long nx, long ny, const double* restrict xbounds, double*
}
}
static
void gridGenYbounds2D(long nx, long ny, const double* restrict ybounds, double* restrict ybounds2D)
void grid_gen_ybounds2D(int nx, int ny, const double *restrict ybounds, double *restrict ybounds2D)
{
long index;
int index;
double minlat, maxlat;
#if defined(_OPENMP)
......@@ -416,7 +414,7 @@ void gridGenYbounds2D(long nx, long ny, const double* restrict ybounds, double*
shared(nx, ny, ybounds, ybounds2D) \
private(minlat, maxlat, index)
#endif
for ( long j = 0; j < ny; ++j )
for ( int j = 0; j < ny; ++j )
{
if ( ybounds[0] > ybounds[1] )
{
......@@ -429,7 +427,7 @@ void gridGenYbounds2D(long nx, long ny, const double* restrict ybounds, double*
minlat = ybounds[2*j ];
}
for ( long i = 0; i < nx; ++i )
for ( int i = 0; i < nx; ++i )
{
index = j*4*nx + 4*i;
ybounds2D[index ] = minlat;
......@@ -945,14 +943,11 @@ int check_range(long n, double *vals, double valid_min, double valid_max)
int gridToCurvilinear(int gridID1, int lbounds)
{
int gridID2;
int gridtype;
size_t gridsize;
long index;
gridtype = gridInqType(gridID1);
gridsize = (size_t) gridInqSize(gridID1);
gridID2 = gridCreate(GRID_CURVILINEAR, (int) gridsize);
int gridtype = gridInqType(gridID1);
size_t gridsize = (size_t) gridInqSize(gridID1);
int gridID2 = gridCreate(GRID_CURVILINEAR, (int) gridsize);
gridDefPrec(gridID2, DATATYPE_FLT32);
switch (gridtype)
......@@ -965,7 +960,6 @@ int gridToCurvilinear(int gridID1, int lbounds)
case GRID_SINUSOIDAL:
{
long i, j;
int nx, ny;
double *xvals = NULL, *yvals = NULL;
double *xvals2D, *yvals2D;
double *xbounds = NULL, *ybounds = NULL;
......@@ -973,8 +967,8 @@ int gridToCurvilinear(int gridID1, int lbounds)
char xunits[CDI_MAX_NAME], yunits[CDI_MAX_NAME];
double xscale = 1, yscale = 1;
nx = gridInqXsize(gridID1);
ny = gridInqYsize(gridID1);
int nx = gridInqXsize(gridID1);
int ny = gridInqYsize(gridID1);
gridInqXunits(gridID1, xunits);
gridInqYunits(gridID1, yunits);
......@@ -1208,8 +1202,8 @@ int gridToCurvilinear(int gridID1, int lbounds)
}
else
{
gridGenXbounds2D(nx, ny, xbounds, xbounds2D);
gridGenYbounds2D(nx, ny, ybounds, ybounds2D);
grid_gen_xbounds2D(nx, ny, xbounds, xbounds2D);
grid_gen_ybounds2D(nx, ny, ybounds, ybounds2D);
}
}
......@@ -1357,8 +1351,8 @@ int gridToUnstructured(int gridID1, int lbounds)
}
else
{
gridGenXbounds2D(nx, ny, xbounds, xbounds2D);
gridGenYbounds2D(nx, ny, ybounds, ybounds2D);
grid_gen_xbounds2D(nx, ny, xbounds, xbounds2D);
grid_gen_ybounds2D(nx, ny, ybounds, ybounds2D);
}
gridDefXbounds(gridID2, xbounds2D);
......
......@@ -22,7 +22,13 @@
void grid_to_radian(const char *units, long nvals, double *restrict values, const char *description);
void grid_to_degree(const char *units, long nvals, double *restrict values, const char *description);
void grid_gen_corners(long n, const double* restrict vals, double* restrict corners);
void grid_gen_corners(int n, const double* restrict vals, double* restrict corners);
void grid_gen_bounds(int n, const double *restrict vals, double *restrict bounds);
void grid_check_lat_borders(int n, double *ybounds);
void grid_gen_xbounds2D(int nx, int ny, const double *restrict xbounds, double *restrict xbounds2D);
void grid_gen_ybounds2D(int nx, int ny, const double *restrict ybounds, double *restrict ybounds2D);
void grid_cell_center_to_bounds_X2D(const char* xunitstr, long xsize, long ysize,
const double* restrict grid_center_lon, double* restrict grid_corner_lon, double dlon);
void grid_cell_center_to_bounds_Y2D(const char* yunitstr, long xsize, long ysize,
......
......@@ -1478,9 +1478,8 @@ int compNlon(int nlat)
static
void gen_grid_lonlat(griddes_t *grid, const char *pline, double inc, double lon1, double lon2, double lat1, double lat2)
{
int nlon, nlat, i;
int gridtype = GRID_LONLAT;
char *endptr;
int lbounds = TRUE;
if ( *pline != 0 )
{
......@@ -1491,25 +1490,112 @@ void gen_grid_lonlat(griddes_t *grid, const char *pline, double inc, double lon1
if ( ! isdigit((int) *pline) && !ispunct((int) *pline) ) return;
endptr = (char *) pline;
char *endptr = (char *) pline;
inc = strtod(pline, &endptr);
if ( *endptr != 0 ) return;
if ( *endptr != 0 )
{
pline = endptr;
if ( *pline == '_' ) pline++;
else return;
if ( *pline == 0 ) return;
if ( *pline == 'c' )
{
gridtype = GRID_CURVILINEAR;
pline++;
if ( *pline == '0' )
{
lbounds = FALSE;
pline++;
}
}
else if ( *pline == 'u' )
{
gridtype = GRID_UNSTRUCTURED;
pline++;
if ( *pline == '0' )
{
lbounds = FALSE;
pline++;
}
}
if ( *pline != 0 ) return;
}
if ( inc < 1e-9 ) inc = 1;
}
grid->type = gridtype;
nlon = (int) ((lon2 - lon1)/inc + 0.5);
nlat = (int) ((lat2 - lat1)/inc + 0.5);
grid->xsize = nlon;
grid->ysize = nlat;
int nlon = (int) ((lon2 - lon1)/inc + 0.5);
int nlat = (int) ((lat2 - lat1)/inc + 0.5);
double *xvals = (double*) Malloc(nlon*sizeof(double));
double *yvals = (double*) Malloc(nlat*sizeof(double));
grid->xvals = (double*) Malloc(grid->xsize*sizeof(double));
grid->yvals = (double*) Malloc(grid->ysize*sizeof(double));
for ( int i = 0; i < nlon; ++i ) xvals[i] = lon1 + inc/2 + i*inc;
for ( int i = 0; i < nlat; ++i ) yvals[i] = lat1 + inc/2 + i*inc;
for ( i = 0; i < nlon; ++i ) grid->xvals[i] = lon1 + inc/2 + i*inc;
for ( i = 0; i < nlat; ++i ) grid->yvals[i] = lat1 + inc/2 + i*inc;
if ( gridtype == GRID_LONLAT )
{
grid->xsize = nlon;
grid->ysize = nlat;
grid->xvals = xvals;
grid->yvals = yvals;
xvals = NULL;
yvals = NULL;
}
else
{
double gridsize = nlon*nlat;
double *xvals2D = (double*) Malloc(gridsize*sizeof(double));
double *yvals2D = (double*) Malloc(gridsize*sizeof(double));
for ( int j = 0; j < nlat; j++ )
for ( int i = 0; i < nlon; i++ )
{
xvals2D[j*nlon+i] = xvals[i];
yvals2D[j*nlon+i] = yvals[j];
}
if ( gridtype == GRID_CURVILINEAR )
{
grid->xsize = nlon;
grid->ysize = nlat;
}
else
{
grid->xsize = gridsize;
grid->ysize = gridsize;
if ( lbounds ) grid->nvertex = 4;
}
grid->xvals = xvals2D;
grid->yvals = yvals2D;
if ( lbounds && nlon > 1 && nlat > 1 )
{
double *xbounds = (double*) Malloc(2*nlon*sizeof(double));
grid_gen_bounds(nlon, xvals, xbounds);
double *ybounds = (double*) Malloc(2*nlat*sizeof(double));
grid_gen_bounds(nlat, yvals, ybounds);
grid_check_lat_borders(2*nlat, ybounds);
double *xbounds2D = (double*) Malloc(4*gridsize*sizeof(double));
double *ybounds2D = (double*) Malloc(4*gridsize*sizeof(double));
grid_gen_xbounds2D(nlon, nlat, xbounds, xbounds2D);
grid_gen_ybounds2D(nlon, nlat, ybounds, ybounds2D);
Free(xbounds);
Free(ybounds);
grid->xbounds = xbounds2D;
grid->ybounds = ybounds2D;
}
}
if ( xvals ) Free(xvals);
if ( yvals ) Free(yvals);
}
......
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