Skip to content
GitLab
Commit 3a21b7a4 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

added grid mask support

parent f27df031
Hide whitespace changes
Inline Side-by-side
src/Consecstat.c
View file @ 3a21b7a4
......@@ -133,7 +133,7 @@ void *Consecstat (void *argument)
int i;
int istreamID, itaxisID, ivlistID, itsID;
int ostreamID, otaxisID, ovlistID, otsID;
int vdate, vtime;
int vdate = 0, vtime = 0;
int histvdate = 0, histvtime = 0;
int recID, nrecs;
int varID, nvars;
......
src/Setgrid.c
View file @ 3a21b7a4
......@@ -20,6 +20,8 @@
Setgrid setgrid Set grid
Setgrid setgridtype Set grid type
Setgrid setgridarea Set grid area
Setgrid setgridmask Set grid mask
*/
#include "cdi.h"
......@@ -32,7 +34,7 @@
void *Setgrid(void *argument)
{
static char func[] = "Setgrid";
int SETGRID, SETGRIDTYPE, SETGRIDAREA;
int SETGRID, SETGRIDTYPE, SETGRIDAREA, SETGRIDMASK;
int operatorID;
int streamID1, streamID2 = CDI_UNDEFID;
int nrecs;
......@@ -41,12 +43,15 @@ void *Setgrid(void *argument)
int taxisID1, taxisID2;
int gridID1, gridID2 = -1;
int ngrids, index;
int gridsize, gridtype = -1;
int gridtype = -1;
int nmiss;
int found;
int areasize = 0;
long i, gridsize;
long areasize = 0;
long masksize = 0;
int lregular = 0;
char *gridname = NULL;
double *gridmask = NULL;
double *areaweight = NULL;
double *array = NULL;
......@@ -55,6 +60,7 @@ void *Setgrid(void *argument)
SETGRID = cdoOperatorAdd("setgrid", 0, 0, "grid description file or name");
SETGRIDTYPE = cdoOperatorAdd("setgridtype", 0, 0, "grid type");
SETGRIDAREA = cdoOperatorAdd("setgridarea", 0, 0, "filename with area weights");
SETGRIDMASK = cdoOperatorAdd("setgridmask", 0, 0, "filename with grid mask");
operatorID = cdoOperatorID();
......@@ -99,7 +105,6 @@ void *Setgrid(void *argument)
if ( cdoVerbose )
{
int i;
double arrmean, arrmin, arrmax;
arrmean = areaweight[0];
......@@ -116,6 +121,33 @@ void *Setgrid(void *argument)
cdoPrint("areaweights: %d %#12.5g%#12.5g%#12.5g", areasize, arrmin, arrmean, arrmax);
}
}
else if ( operatorID == SETGRIDMASK )
{
int streamID, vlistID, gridID;
char *maskfile;
double missval;
maskfile = operatorArgv()[0];
streamID = streamOpenRead(maskfile);
if ( streamID < 0 ) cdiError(streamID, "Open failed on %s", maskfile);
vlistID = streamInqVlist(streamID);
nrecs = streamInqTimestep(streamID, 0);
streamInqRecord(streamID, &varID, &levelID);
missval = vlistInqVarMissval(vlistID, varID);
gridID = vlistInqVarGrid(vlistID, varID);
masksize = gridInqSize(gridID);
gridmask = (double *) malloc(masksize*sizeof(double));
streamReadRecord(streamID, gridmask, &nmiss);
streamClose(streamID);
for ( i = 0; i < masksize; i++ )
if ( DBL_IS_EQUAL(gridmask[i], missval) ) gridmask[i] = 0;
}
streamID1 = streamOpenRead(cdoStreamName(0));
if ( streamID1 < 0 ) cdiError(streamID1, "Open failed on %s", cdoStreamName(0));
......@@ -187,6 +219,32 @@ void *Setgrid(void *argument)
}
}
}
else if ( operatorID == SETGRIDMASK )
{
ngrids = vlistNgrids(vlistID1);
for ( index = 0; index < ngrids; index++ )
{
gridID1 = vlistGrid(vlistID1, index);
gridtype = gridInqType(gridID1);
gridsize = gridInqSize(gridID1);
if ( gridsize == masksize )
{
int *mask;
mask = (int *) malloc(masksize*sizeof(int));
for ( i = 0; i < masksize; i++ )
{
if ( gridmask[i] < 0 || gridmask[i] > 255 )
mask[i] = 0;
else
mask[i] = NINT(gridmask[i]);
}
gridID2 = gridDuplicate(gridID1);
gridDefMask(gridID2, mask);
vlistChangeGridIndex(vlistID2, index, gridID2);
free(mask);
}
}
}
streamDefVlist(streamID2, vlistID2);
//vlistPrint(vlistID2);
......@@ -230,6 +288,7 @@ void *Setgrid(void *argument)
streamClose(streamID2);
streamClose(streamID1);
if ( gridmask ) free(gridmask);
if ( areaweight ) free(areaweight);
if ( array ) free(array);
......
src/Writegrid.c
View file @ 3a21b7a4
......@@ -31,118 +31,48 @@
void *Writegrid(void *argument)
{
static char func[] = "Writegrid";
int WRITEGRID, GRIDAREA;
int streamID;
int vlistID1;
int vlistID;
int gridID;
int operatorID;
int i;
int gridtype, gridsize;
int varID, levelID;
int nrecs;
int nmiss;
int *imask = NULL;
double missval;
double *array = NULL;
long gridsize, i;
int gridtype;
int *mask = NULL;
cdoInitialize(argument);
WRITEGRID = cdoOperatorAdd("writegrid", 0, 0, NULL);
GRIDAREA = cdoOperatorAdd("gridarea", 0, 0, NULL);
operatorID = cdoOperatorID();
streamID = streamOpenRead(cdoStreamName(0));
if ( streamID < 0 ) cdiError(streamID, "Open failed on %s", cdoStreamName(0));
vlistID1 = streamInqVlist(streamID);
if ( operatorID == WRITEGRID )
{
nrecs = streamInqTimestep(streamID, 0);
streamInqRecord(streamID, &varID, &levelID);
vlistID = streamInqVlist(streamID);
gridID = vlistGrid(vlistID, 0);
gridID = vlistInqVarGrid(vlistID1, varID);
gridtype = gridInqType(gridID);
gridsize = gridInqSize(gridID);
gridtype = gridInqType(gridID);
gridsize = gridInqSize(gridID);
if ( gridtype == GRID_GME ) gridID = gridToCell(gridID);
if ( gridtype == GRID_GME ) gridID = gridToCell(gridID);
if ( gridtype != GRID_CURVILINEAR && gridtype != GRID_CELL )
gridID = gridToCurvilinear(gridID);
if ( gridtype != GRID_CURVILINEAR && gridtype != GRID_CELL )
gridID = gridToCurvilinear(gridID);
if ( gridInqXbounds(gridID, NULL) == 0 || gridInqYbounds(gridID, NULL) == 0 )
cdoAbort("Grid corner missing!");
if ( gridInqXbounds(gridID, NULL) == 0 || gridInqYbounds(gridID, NULL) == 0 )
cdoAbort("Grid corner missing!");
array = (double *) malloc(gridsize*sizeof(double));
imask = (int *) malloc(gridsize*sizeof(int));
streamReadRecord(streamID, array, &nmiss);
mask = (int *) malloc(gridsize*sizeof(int));
missval = vlistInqVarMissval(vlistID1, varID);
for ( i = 0; i < gridsize; i++ )
{
if ( DBL_IS_EQUAL(array[i], missval) )
imask[i] = 0;
else
imask[i] = 1;
}
writeNCgrid(cdoStreamName(1), gridID, imask);
if ( gridInqMask(gridID, NULL) )
{
gridInqMask(gridID, mask);
}
else if ( operatorID == GRIDAREA )
else
{
double *area, *xvals, *yvals;
int zaxisID, vlistID2, streamID2;
int nlon, nlat;
gridID = vlistGrid(vlistID1, 0);
gridtype = gridInqType(gridID);
gridsize = gridInqSize(gridID);
if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_LONLAT )
{
area = (double *) malloc(gridsize*sizeof(double));
xvals = (double *) malloc(gridsize*sizeof(double));
yvals = (double *) malloc(gridsize*sizeof(double));
zaxisID = zaxisCreate(ZAXIS_SURFACE, 1);
vlistID2 = vlistCreate();
varID = vlistDefVar(vlistID2, gridID, zaxisID, TIME_CONSTANT);
/* vlistDefVarCode(vlistID2, varID, 11); */
vlistDefVarName(vlistID2, varID, "area");
vlistDefVarLongname(vlistID2, varID, "area");
vlistDefVarUnits(vlistID2, varID, "m**2");
gridID = gridToCurvilinear(gridID);
nlon = gridInqXsize(gridID);
nlat = gridInqYsize(gridID);
gridInqXvals(gridID, xvals);
gridInqYvals(gridID, yvals);
streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
streamDefVlist(streamID2, vlistID2);
streamDefRecord(streamID2, varID, 0);
nmiss = 0;
memcpy(area, xvals, gridsize*sizeof(double));
streamWriteRecord(streamID2, area, nmiss);
streamClose(streamID2);
free(area);
free(xvals);
free(yvals);
}
else
cdoAbort("Unsupported grid type: %s", gridNamePtr(gridtype));
for ( i = 0; i < gridsize; i++ ) mask[i] = 1;
}
writeNCgrid(cdoStreamName(1), gridID, mask);
streamClose(streamID);
if (array) free(array);
if (imask) free(imask);
if ( mask ) free(mask);
cdoFinish();
......
src/grid.c
View file @ 3a21b7a4
......@@ -560,6 +560,21 @@ int gridToRegular(int gridID1)
return (gridID2);
}
static
void gridCopyMask(int gridID1, int gridID2, long gridsize)
{
static char func[] = "gridCopyMask";
if ( gridInqMask(gridID1, NULL) )
{
int *mask;
mask = (int *) malloc(gridsize*sizeof(int));
gridInqMask(gridID1, mask);
gridDefMask(gridID2, mask);
free(mask);
}
}
int gridToCurvilinear(int gridID1)
{
......@@ -803,6 +818,8 @@ int gridToCurvilinear(int gridID1)
}
}
gridCopyMask(gridID1, gridID2, gridsize);
break;
}
default:
......@@ -932,6 +949,8 @@ int gridToCell(int gridID1)
free(ybounds2D);
}
gridCopyMask(gridID1, gridID2, gridsize);
break;
}
case GRID_GME:
......@@ -991,6 +1010,8 @@ int gridToCell(int gridID1)
free (xbounds);
free (ybounds);
gridCopyMask(gridID1, gridID2, gridsize);
break;
}
default:
......@@ -1003,7 +1024,6 @@ int gridToCell(int gridID1)
return (gridID2);
}
static
double areas(struct cart *dv1, struct cart *dv2, struct cart *dv3)
{
......@@ -1090,7 +1110,6 @@ double areas(struct cart *dv1, struct cart *dv2, struct cart *dv3)
return areas;
}
static
double cell_area(long i, long nv, double *grid_center_lon, double *grid_center_lat,
double *grid_corner_lon, double *grid_corner_lat)
......
src/griddes.c
View file @ 3a21b7a4
......@@ -58,6 +58,7 @@ void extClose(int fileID);
void gridInit(grid_t *grid)
{
grid->mask = NULL;
grid->xvals = NULL;
grid->yvals = NULL;
grid->xbounds = NULL;
......@@ -270,6 +271,13 @@ int gridDefine(grid_t grid)
gridDefXpole(gridID, grid.xpole);
gridDefYpole(gridID, grid.ypole);
}
if ( grid.mask )
{
gridDefMask(gridID, grid.mask);
free(grid.mask);
}
break;
}
case GRID_CURVILINEAR:
......@@ -323,6 +331,12 @@ int gridDefine(grid_t grid)
free(grid.ybounds);
}
if ( grid.mask )
{
gridDefMask(gridID, grid.mask);
free(grid.mask);
}
break;
}
case GRID_LCC:
......@@ -350,6 +364,12 @@ int gridDefine(grid_t grid)
gridDefLCC(gridID, grid.originLon, grid.originLat, grid.lonParY,
grid.lat1, grid.lat2, grid.xinc, grid.yinc, grid.projflag, grid.scanflag);
if ( grid.mask )
{
gridDefMask(gridID, grid.mask);
free(grid.mask);
}
break;
}
case GRID_LCC2:
......@@ -399,6 +419,12 @@ int gridDefine(grid_t grid)
gridDefLcc2(gridID, grid.a, grid.lon_0, grid.lat_0, grid.lat_1, grid.lat_2);
if ( grid.mask )
{
gridDefMask(gridID, grid.mask);
free(grid.mask);
}
break;
}
case GRID_SPECTRAL:
......@@ -415,13 +441,13 @@ int gridDefine(grid_t grid)
gridDefTrunc(gridID, grid.ntr);
gridDefComplexPacking(gridID, grid.lcomplex);
break;
}
case GRID_GME:
{
if ( grid.nd == 0 ) Error(func, "nd undefined!");
if ( grid.ni == 0 ) Error(func, "ni undefined!");
if ( grid.nd == 0 ) Error(func, "nd undefined!");
if ( grid.ni == 0 ) Error(func, "ni undefined!");
if ( grid.size == 0 ) Error(func, "size undefined!");
gridID = gridCreate(grid.type, grid.size);
......@@ -433,6 +459,12 @@ int gridDefine(grid_t grid)
gridDefGMEni2(gridID, grid.ni2);
gridDefGMEni3(gridID, grid.ni3);
if ( grid.mask )
{
gridDefMask(gridID, grid.mask);
free(grid.mask);
}
break;
}
default:
......@@ -454,7 +486,8 @@ int gridDefine(grid_t grid)
return (gridID);
}
static char *skipSeparator(char *pline)
static
char *skipSeparator(char *pline)
{
while ( isspace((int) *pline) ) pline++;
if ( *pline == '=' || *pline == ':' ) pline++;
......@@ -463,6 +496,7 @@ static char *skipSeparator(char *pline)
return (pline);
}
void fnmexp2(char *out, char *in1, const char *in2)
{
const char *pos, *ch;
......@@ -893,6 +927,40 @@ int gridFromFile(FILE *gfp, const char *dname)
grid.xvals[i] = flon;
}
}
else if ( cmpstr(pline, "mask", len) == 0 )
{
int i = 0;
long lval;
char *endptr;
size = grid.size;
if ( size > 0 )
{
pline = skipSeparator(pline + len);
grid.mask = (int *) malloc(size*sizeof(int));
for ( i = 0; i < size; i++ )
{
endptr = pline;
lval = strtol(pline, &endptr, 10);
if ( pline == endptr )
{
if ( ! readline(gfp, line, MAX_LINE_LEN) )
{
Warning(func, "Incomplete command: >mask<");
break;
}
pline = line;
lval = strtol(pline, &endptr, 10);
}
grid.mask[i] = (int)lval;
pline = endptr;
}
}
else
Warning(func, "gridsize undefined!");
}
else if ( cmpstr(pline, "xvals", len) == 0 )
{
int i = 0;
......
src/griddes.h
View file @ 3a21b7a4
......@@ -2,6 +2,7 @@
#define _GRIDDES_H
typedef struct {
int *mask;
double *xvals;
double *yvals;
double *xbounds;
......
src/griddes_nc.c
View file @ 3a21b7a4
......@@ -57,6 +57,7 @@ int gridFromNCfile(const char *gridfile)
int nc_gridclon_id; /* netCDF grid corner lon var id */
int nc_gridlat_id; /* netCDF grid center lat var id */
int nc_gridlon_id; /* netCDF grid center lon var id */
int nc_gridmask_id; /* netCDF grid mask id */
nc_type xtype;
size_t attlen;
......@@ -86,7 +87,6 @@ int gridFromNCfile(const char *gridfile)
nc_inq_varid(nc_file_id, "grid_corner_lat", &nc_gridclat_id) != NC_NOERR ||
nc_inq_varid(nc_file_id, "grid_corner_lon", &nc_gridclon_id) != NC_NOERR ) return (gridID);
nce(nc_inq_varid(nc_file_id, "grid_dims", &nc_griddims_id));
nce(nc_get_var_int(nc_file_id, nc_griddims_id, grid_dims));
if ( grid_rank == 1 )
......@@ -112,11 +112,6 @@ int gridFromNCfile(const char *gridfile)
grid.xbounds = (double *) malloc(grid.nvertex*grid.size*sizeof(double));
grid.ybounds = (double *) malloc(grid.nvertex*grid.size*sizeof(double));
nce(nc_inq_varid(nc_file_id, "grid_center_lat", &nc_gridlat_id));
nce(nc_inq_varid(nc_file_id, "grid_center_lon", &nc_gridlon_id));
nce(nc_inq_varid(nc_file_id, "grid_corner_lat", &nc_gridclat_id));
nce(nc_inq_varid(nc_file_id, "grid_corner_lon", &nc_gridclon_id));
nce(nc_inq_vartype(nc_file_id, nc_gridlat_id, &xtype));
if ( xtype == NC_FLOAT ) grid.prec = DATATYPE_FLT32;
else grid.prec = DATATYPE_FLT64;
......@@ -133,6 +128,21 @@ int gridFromNCfile(const char *gridfile)
nce(nc_get_att_text(nc_file_id, nc_gridlat_id, "units", grid.yunits));
grid.yunits[attlen] = 0;
if ( nc_inq_varid(nc_file_id, "grid_imask", &nc_gridmask_id) == NC_NOERR )
{
int i;
grid.mask = (int *) malloc(grid.size*sizeof(int));
nce(nc_get_var_int(nc_file_id, nc_gridmask_id, grid.mask));
for ( i = 0; i < grid.size; ++i )
if ( grid.mask[i] != 1 ) break;
if ( i == grid.size )
{
free(grid.mask);
grid.mask = NULL;
}
}
gridID = gridDefine(grid);
}
......
src/modules.c
View file @ 3a21b7a4
......@@ -345,7 +345,7 @@ void *Wct(void *argument);
#define SetOperators {"setpartab", "setpartabv", "setcode", "setname", "setlevel", "setltype", "settabnum"}
#define SetboxOperators {"setclonlatbox", "setcindexbox"}
#define SetgattOperators {"setgatt", "setgatts"}
#define SetgridOperators {"setgrid", "setgridtype", "setgridarea"}
#define SetgridOperators {"setgrid", "setgridtype", "setgridarea", "setgridmask"}
#define SethaloOperators {"sethalo", "tpnhalo"}
#define SetmissOperators {"setmissval", "setctomiss", "setmisstoc", "setrtomiss", "setvrange"}
#define SetrcanameOperators {"setrcaname"}
......
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