Commit 0fa7faa6 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

selgridcell: read indices from mask.

parent e305b290
......@@ -3,6 +3,10 @@
* Using CDI library version 1.9.0
* Version 1.9.0 release
2017-06-04 Uwe Schulzweida
* selgridcell: read indices from mask file
2017-06-02 Uwe Schulzweida
* New operator ensrange: Ensemble range (ensmax-ensmin)
......
......@@ -65,11 +65,9 @@ void sel_index(double *array1, double *array2, int nind, int *indarr)
void *Selgridcell(void *argument)
{
int nrecs;
int varID, levelID;
int varID;
int gridID1 = -1, gridID2;
int index, gridtype = -1;
int nmiss;
double missval;
typedef struct {
int gridID1, gridID2;
} sindex_t;
......@@ -77,17 +75,42 @@ void *Selgridcell(void *argument)
cdoInitialize(argument);
int SELGRIDCELL = cdoOperatorAdd("selgridcell", 0, 0, "grid cell indices (1-N)");
cdoOperatorAdd("selgridcell", 0, 0, "grid cell indices (1-N)");
int DELGRIDCELL = cdoOperatorAdd("delgridcell", 0, 0, "grid cell indices (1-N)");
UNUSED(SELGRIDCELL);
operatorInputArg(cdoOperatorEnter(0));
int operatorID = cdoOperatorID();
int nind = args2int_lista(operatorArgc(), operatorArgv(), ilista);
int *indarr = (int*) lista_dataptr(ilista);
int nind;
int *indarr = NULL;
if ( operatorArgc() == 1 && fileExists(operatorArgv()[0]) )
{
bool *cdo_read_mask(const char *maskfile, int *n);
int n = 0;
bool *mask = cdo_read_mask(operatorArgv()[0], &n);
nind = 0;
for ( int i = 0; i < n; ++i ) if ( mask[i] ) nind++;
if ( nind == 0 ) cdoAbort("Mask is empty!");
else
{
indarr = (int*) Malloc(nind*sizeof(double));
nind = 0;
for ( int i = 0; i < n; ++i ) if ( mask[i] ) indarr[nind++] = i;
}
if ( mask ) Free(mask);
}
else
{
nind = args2int_lista(operatorArgc(), operatorArgv(), ilista);
indarr = (int*) lista_dataptr(ilista);
if ( cdoVerbose )
for ( int i = 0; i < nind; i++ )
cdoPrint("int %d = %d", i+1, indarr[i]);
for ( int i = 0; i < nind; i++ ) indarr[i] -= 1;
}
int indmin = indarr[0];
int indmax = indarr[0];
......@@ -97,14 +120,7 @@ void *Selgridcell(void *argument)
if ( indmin > indarr[i] ) indmin = indarr[i];
}
if ( cdoVerbose )
for ( int i = 0; i < nind; i++ )
cdoPrint("int %d = %d", i+1, indarr[i]);
if ( indmin < 1 ) cdoAbort("Index < 1 not allowed!");
for ( int i = 0; i < nind; i++ ) indarr[i] -= 1;
if ( indmin < 0 ) cdoAbort("Index < 1 not allowed!");
int streamID1 = streamOpenRead(cdoStreamName(0));
......@@ -136,6 +152,7 @@ void *Selgridcell(void *argument)
if ( cellidx[i] == 1 ) cellidx[j++] = i;
if ( j != ncells ) cdoAbort("Internal error; number of cells differ");
}
if ( ncells == 0 ) cdoAbort("Mask is empty!");
for ( index = 0; index < ngrids; index++ )
{
......@@ -144,7 +161,7 @@ void *Selgridcell(void *argument)
int gridsize = gridInqSize(gridID1);
if ( gridsize == 1 ) continue;
if ( indmax > gridsize )
if ( indmax >= gridsize )
{
cdoWarning("Max grid index is greater than grid size, skipped grid %d!", index+1);
continue;
......@@ -195,6 +212,7 @@ void *Selgridcell(void *argument)
for ( int recID = 0; recID < nrecs; recID++ )
{
int nmiss, varID, levelID;
streamInqRecord(streamID1, &varID, &levelID);
streamReadRecord(streamID1, array1, &nmiss);
......@@ -216,7 +234,7 @@ void *Selgridcell(void *argument)
if ( nmiss )
{
nmiss = 0;
missval = vlistInqVarMissval(vlistID2, varID);
double missval = vlistInqVarMissval(vlistID2, varID);
for ( int i = 0; i < gridsize2; i++ )
if ( DBL_IS_EQUAL(array2[i], missval) ) nmiss++;
}
......
#include <cdi.h>
#include "cdo_int.h"
bool *cdo_read_timestepmask(const char *maskfile, int *n)
{
*n = 0;
......@@ -58,3 +59,43 @@ bool *cdo_read_timestepmask(const char *maskfile, int *n)
return imask;
}
bool *cdo_read_mask(const char *maskfile, int *n)
{
*n = 0;
int streamID = streamOpenRead(maskfile);
if ( streamID == CDI_UNDEFID ) cdoAbort("Open failed on %s!", maskfile);
int vlistID = streamInqVlist(streamID);
int nvars = vlistNvars(vlistID);
if ( nvars > 1 ) cdoAbort("Mask %s contains more than one variable!", maskfile);
int gridsize = gridInqSize(vlistInqVarGrid(vlistID, 0));
int nlev = zaxisInqSize(vlistInqVarZaxis(vlistID, 0));
if ( nlev > 1 ) cdoAbort("Mask %s has more than one level!", maskfile);
*n = gridsize;
bool *imask = (bool*) Malloc(gridsize*sizeof(bool));
double *dmask = (double*) Malloc(gridsize*sizeof(double));
int nrecs = streamInqTimestep(streamID, 0);
if ( nrecs != 1 ) cdoAbort("Internal error; unexprected number of records!");
int varID, levelID;
int nmiss;
streamInqRecord(streamID, &varID, &levelID);
streamReadRecord(streamID, dmask, &nmiss);
for ( int i = 0; i < gridsize; ++i )
imask[i] = IS_NOT_EQUAL(dmask[i], 0);
Free(dmask);
streamClose(streamID);
return imask;
}
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