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

added function gridmask

parent 81e339b7
......@@ -61,7 +61,7 @@ dimensions: ...
variables: ...
// global attributes:
:myatt = "myattcontens" ;
:myatt = "myattcontents" ;
}
@EndVerbatim
@EndExample
......@@ -2,7 +2,7 @@
This file is part of CDO. CDO is a collection of Operators to
manipulate and analyse Climate model Data.
Copyright (C) 2003-2009 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
Copyright (C) 2003-2010 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
See COPYING file for copying and redistribution conditions.
This program is free software; you can redistribute it and/or modify
......@@ -20,6 +20,7 @@
Gridcell gridarea Grid cell area in m^2
Gridcell gridweights Grid cell weights
Gridcell gridmask Grid mask
*/
......@@ -32,19 +33,18 @@
void *Gridcell(void *argument)
{
static char func[] = "Gridcell";
int GRIDAREA, GRIDWGTS;
int GRIDAREA, GRIDWGTS, GRIDMASK;
int operatorID;
int streamID1, streamID2;
int vlistID1, vlistID2;
int gridID, zaxisID;
int gridsize, gridtype;
int gridtype;
int status;
int ngrids;
int tsID, varID, levelID, taxisID;
double *grid_area = NULL;
double *grid_wgts = NULL;
double *pdata;
long i, gridsize;
char *envstr;
double *array = NULL;
double EarthRadius = 6371000; /* default radius of the earth in m */
double PlanetRadius = EarthRadius;
......@@ -52,6 +52,7 @@ void *Gridcell(void *argument)
GRIDAREA = cdoOperatorAdd("gridarea", 0, 0, NULL);
GRIDWGTS = cdoOperatorAdd("gridweights", 0, 0, NULL);
GRIDMASK = cdoOperatorAdd("gridmask", 0, 0, NULL);
operatorID = cdoOperatorID();
......@@ -96,7 +97,14 @@ void *Gridcell(void *argument)
vlistDefVarUnits(vlistID2, varID, "m2");
}
else if ( operatorID == GRIDWGTS )
vlistDefVarName(vlistID2, varID, "cell_weights");
{
vlistDefVarName(vlistID2, varID, "cell_weights");
}
else if ( operatorID == GRIDMASK )
{
vlistDefVarName(vlistID2, varID, "grid_mask");
vlistDefVarDatatype(vlistID2, varID, DATATYPE_UINT8);
}
taxisID = taxisCreate(TAXIS_ABSOLUTE);
vlistDefTaxis(vlistID2, taxisID);
......@@ -108,8 +116,7 @@ void *Gridcell(void *argument)
gridsize = gridInqSize(gridID);
grid_area = (double *) malloc(gridsize*sizeof(double));
grid_wgts = (double *) malloc(gridsize*sizeof(double));
array = (double *) malloc(gridsize*sizeof(double));
tsID = 0;
......@@ -133,20 +140,18 @@ void *Gridcell(void *argument)
if ( gridHasArea(gridID) )
{
if ( cdoVerbose ) cdoPrint("Using existing grid cell area!");
gridInqArea(gridID, grid_area);
gridInqArea(gridID, array);
}
else
{
int i;
status = gridGenArea(gridID, grid_area);
status = gridGenArea(gridID, array);
if ( status == 1 )
cdoAbort("Grid corner missing!");
else if ( status == 2 )
cdoAbort("Can't compute grid cell areas for this grid!");
for ( i = 0; i < gridsize; ++i )
grid_area[i] *= PlanetRadius*PlanetRadius;
array[i] *= PlanetRadius*PlanetRadius;
}
}
else
......@@ -157,26 +162,37 @@ void *Gridcell(void *argument)
else
cdoAbort("Unsupported grid type: %s", gridNamePtr(gridtype));
}
pdata = grid_area;
}
else
else if ( operatorID == GRIDWGTS )
{
status = gridWeights(gridID, grid_wgts);
status = gridWeights(gridID, array);
if ( status != 0 )
cdoWarning("Using constant grid cell area weights!");
}
else if ( operatorID == GRIDMASK )
{
int *mask;
mask = (int *) malloc(gridsize*sizeof(int));
if ( gridInqMask(gridID, NULL) )
{
gridInqMask(gridID, mask);
}
else
{
for ( i = 0; i < gridsize; ++i ) mask[i] = 1;
}
pdata = grid_wgts;
for ( i = 0; i < gridsize; ++i ) array[i] = mask[i];
free(mask);
}
streamWriteRecord(streamID2, pdata, 0);
streamWriteRecord(streamID2, array, 0);
streamClose(streamID2);
streamClose(streamID1);
if ( grid_area ) free(grid_area);
if ( grid_wgts ) free(grid_wgts);
if ( array ) free(array);
cdoFinish();
......
......@@ -650,10 +650,10 @@ void *Outputgmt(void *argument)
nlat = gridInqYsize(gridID);
nlev = zaxisInqSize(zaxisID);
if ( gridInqMask(gridID, NULL) )
if ( gridInqMaskGME(gridID, NULL) )
{
grid_mask = (int *) malloc(gridsize*sizeof(int));
gridInqMask(gridID, grid_mask);
gridInqMaskGME(gridID, grid_mask);
}
if ( gridInqType(gridID) != GRID_CELL )
......
......@@ -534,7 +534,7 @@ void *Remap(void *argument)
remaps[0].grid.grid2_vgpm = (int *) realloc(remaps[0].grid.grid2_vgpm,
gridInqSize(gridID2)*sizeof(int));
gridID2_gme = gridToCell(gridID2);
gridInqMask(gridID2_gme, remaps[0].grid.grid2_vgpm);
gridInqMaskGME(gridID2_gme, remaps[0].grid.grid2_vgpm);
for ( i = 0; i < gridsize2; ++i )
if ( remaps[0].grid.grid2_vgpm[i] ) isize++;
gridsize2 = isize;
......
......@@ -490,6 +490,9 @@ char *gridNamePtr(int gridtype);
void gridCompress(int gridID);
void gridDefMaskGME(int gridID, const int *mask_vec);
int gridInqMaskGME(int gridID, int *mask_vec);
void gridDefMask(int gridID, const int *mask_vec);
int gridInqMask(int gridID, int *mask_vec);
......
......@@ -975,7 +975,7 @@ int gridToCell(int gridID1)
gridDefXvals(gridID2, xvals);
gridDefYvals(gridID2, yvals);
gridDefMask(gridID2, imask);
gridDefMaskGME(gridID2, imask);
gridDefNvertex(gridID2, nv);
......@@ -1201,7 +1201,7 @@ int gridGenArea(int gridID, double *area)
lgriddestroy = TRUE;
gridID = gridToCell(gridID);
grid_mask = (int *) malloc(gridsize*sizeof(int));
gridInqMask(gridID, grid_mask);
gridInqMaskGME(gridID, grid_mask);
}
else
{
......@@ -1294,7 +1294,7 @@ int gridGenWeights(int gridID, double *grid_area, double *grid_wgts)
{
gridID = gridToCell(gridID);
grid_mask = (int *) malloc(gridsize*sizeof(int));
gridInqMask(gridID, grid_mask);
gridInqMaskGME(gridID, grid_mask);
}
total_area = 0;
......
......@@ -281,7 +281,7 @@ void *Wct(void *argument);
#define GradsdesOperators {"gradsdes1", "gradsdes2", "dumpmap"}
#define GridboxstatOperators {"gridboxmin", "gridboxmax", "gridboxsum", "gridboxmean", "gridboxavg", \
"gridboxvar", "gridboxstd"}
#define GridcellOperators {"gridarea", "gridweights"}
#define GridcellOperators {"gridarea", "gridweights", "gridmask"}
#define HarmonicOperators {"harmonic"}
#define HistogramOperators {"histcount", "histsum", "histmean", "histfreq"}
#define ImportamsrOperators {"import_amsr"}
......
......@@ -914,7 +914,7 @@ void remapGridInit(int map_type, int lextrapolate, int gridID1, int gridID2, rem
#endif
for ( i = 0; i < grid1_size; i++ ) rg->grid1_mask[i] = TRUE;
if ( gridInqType(rg->gridID1) == GRID_GME ) gridInqMask(gridID1_gme, rg->grid1_vgpm);
if ( gridInqType(rg->gridID1) == GRID_GME ) gridInqMaskGME(gridID1_gme, rg->grid1_vgpm);
/* Convert lat/lon units if required */
......@@ -975,7 +975,7 @@ void remapGridInit(int map_type, int lextrapolate, int gridID1, int gridID2, rem
#endif
for ( i = 0; i < grid2_size; i++ ) rg->grid2_mask[i] = TRUE;
if ( gridInqType(rg->gridID2) == GRID_GME ) gridInqMask(gridID2_gme, rg->grid2_vgpm);
if ( gridInqType(rg->gridID2) == GRID_GME ) gridInqMaskGME(gridID2_gme, rg->grid2_vgpm);
/* Convert lat/lon units if required */
......@@ -6761,7 +6761,7 @@ void read_remap_scrip(const char *interp_file, int gridID1, int gridID2, int *ma
remapGridRealloc(rv->map_type, rg);
if ( gridInqType(gridID1) == GRID_GME ) gridInqMask(gridID1_gme_c, rg->grid1_vgpm);
if ( gridInqType(gridID1) == GRID_GME ) gridInqMaskGME(gridID1_gme_c, rg->grid1_vgpm);
rv->pinit = TRUE;
for ( i = 0; i < 4; i++ ) rv->wts[i] = NULL;
......
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