Skip to content

Commit

Permalink
apps/radar/Ecco - adding terrain ht and water layer
Browse files Browse the repository at this point in the history
  • Loading branch information
@mike-dixon
mike-dixon committed Apr 3, 2024
1 parent e324b1b commit 000ee40
Show file tree
Hide file tree
Showing 6 changed files with 164 additions and 80 deletions.
170 changes: 115 additions & 55 deletions codebase/apps/radar/src/Ecco/Ecco.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@
#include "TerrainHt.hh"
using namespace std;

const fl32 Ecco::_missing = -9999.0;
const fl32 Ecco::_missingFl32 = -9999.0;

// Constructor

Expand Down Expand Up @@ -349,78 +349,118 @@ int Ecco::_addTerrainHtField()
// get geometry from DBZ field

MdvxField *dbzField = _inMdvx.getField(_params.dbz_field_name);

// create terrain ht field - this is 2D

Mdvx::field_header_t fhdrTerrain = dbzField->getFieldHeader();
fhdrTerrain.nz = 1;
fhdrTerrain.vlevel_type = Mdvx::VERT_TYPE_SURFACE;
size_t planeSize32 = fhdrTerrain.nx * fhdrTerrain.ny * sizeof(fl32);
fhdrTerrain.volume_size = planeSize32;
fhdrTerrain.encoding_type = Mdvx::ENCODING_FLOAT32;
fhdrTerrain.data_element_nbytes = 4;
fhdrTerrain.missing_data_value = _missing;
fhdrTerrain.bad_data_value = _missing;
fhdrTerrain.scale = 1.0;
fhdrTerrain.bias = 0.0;

Mdvx::vlevel_header_t vhdr2d;
MEM_zero(vhdr2d);
vhdr2d.level[0] = 0;
vhdr2d.type[0] = Mdvx::VERT_TYPE_SURFACE;

MdvxField::setFieldName("TerrainHt", fhdrTerrain);
MdvxField::setFieldNameLong("Terrain height MSL", fhdrTerrain);
MdvxField::setUnits("m", fhdrTerrain);
MdvxField *terrainHtField =
new MdvxField(fhdrTerrain, vhdr2d, NULL, false, false, false);
Mdvx::field_header_t fhdrDbz = dbzField->getFieldHeader();

// create projection object

MdvxProj proj(fhdrTerrain);
MdvxProj proj(fhdrDbz);

// load up height data
// load up height and water data

float *ht = new float[fhdrTerrain.nx * fhdrTerrain.ny];
float *ht = new float[fhdrDbz.nx * fhdrDbz.ny];
ui08 *water = new ui08[fhdrDbz.nx * fhdrDbz.ny];

int ii = 0;
double yy = fhdrTerrain.grid_miny;
for (int iy = 0; iy < fhdrTerrain.ny; iy++, yy += fhdrTerrain.grid_dy) {
double xx = fhdrTerrain.grid_minx;
for (int ix = 0; ix < fhdrTerrain.nx; ix++, xx += fhdrTerrain.grid_dx, ii++) {
int index = 0;
double yy = fhdrDbz.grid_miny;
for (int iy = 0; iy < fhdrDbz.ny; iy++, yy += fhdrDbz.grid_dy) {
double xx = fhdrDbz.grid_minx;
for (int ix = 0; ix < fhdrDbz.nx; ix++, xx += fhdrDbz.grid_dx, index++) {

// get lat/lon
double lat, lon;
proj.xy2latlon(xx, yy, lat, lon);

// get height

double terrainHtM = 0.0;
bool isWater = true;
if (_terrainHt->getHt(lat, lon, terrainHtM, isWater) == 0) {
ht[ii] = terrainHtM;
ht[index] = terrainHtM;
if (isWater) {
water[index] = 1;
} else {
water[index] = 0;
}
} else {
ht[ii] = 0.0;
ht[index] = 0.0;
water[index] = 2;
}

if (_params.debug >= Params::DEBUG_EXTRA) {
cerr << "xx, yy, lat, lon, ht: "
<< xx << ", "
<< yy << ", "
<< lat << ", "
<< lon << ", "
<< terrainHtM << endl;
cerr << "xx, yy, lat, lon, ht, water: "
<< xx << ", " << yy << ", "
<< lat << ", " << lon << ", "
<< terrainHtM << ", " << (isWater?"Y":"N") << endl;
}

} // ix
} // iy

terrainHtField->setVolData(ht, planeSize32, Mdvx::ENCODING_FLOAT32);
// create terrain ht 2D field

Mdvx::field_header_t fhdrTerrain = dbzField->getFieldHeader();
fhdrTerrain.nz = 1;
fhdrTerrain.vlevel_type = Mdvx::VERT_TYPE_SURFACE;
size_t planeSize32 = fhdrTerrain.nx * fhdrTerrain.ny * sizeof(fl32);
fhdrTerrain.volume_size = planeSize32;
fhdrTerrain.encoding_type = Mdvx::ENCODING_FLOAT32;
fhdrTerrain.data_element_nbytes = 4;
fhdrTerrain.missing_data_value = _missingFl32;
fhdrTerrain.bad_data_value = _missingFl32;
fhdrTerrain.scale = 1.0;
fhdrTerrain.bias = 0.0;

Mdvx::vlevel_header_t vhdrTerrain;
MEM_zero(vhdrTerrain);
vhdrTerrain.level[0] = 0;
vhdrTerrain.type[0] = Mdvx::VERT_TYPE_SURFACE;

MdvxField::setFieldName("TerrainHt", fhdrTerrain);
MdvxField::setFieldNameLong("Terrain height MSL", fhdrTerrain);
MdvxField::setUnits("m", fhdrTerrain);
MdvxField *terrainHtField =
new MdvxField(fhdrTerrain, vhdrTerrain, ht, false, false, false);
delete[] ht;

// add field to input object
// add ht field to input object

_inMdvx.addField(terrainHtField);

// return now if no water layer to be added

if (!_params.add_water_layer) {
delete[] water;
return 0;
}

// create water layer field - this is 2D

Mdvx::field_header_t fhdrWater = dbzField->getFieldHeader();
fhdrWater.nz = 1;
fhdrWater.vlevel_type = Mdvx::VERT_TYPE_SURFACE;
size_t planeSize08 = fhdrWater.nx * fhdrWater.ny * sizeof(ui08);
fhdrWater.volume_size = planeSize08;
fhdrWater.encoding_type = Mdvx::ENCODING_INT8;
fhdrWater.data_element_nbytes = 1;
fhdrWater.missing_data_value = 2;
fhdrWater.bad_data_value = 2;
fhdrWater.scale = 1.0;
fhdrWater.bias = 0.0;

Mdvx::vlevel_header_t vhdrWater;
MEM_zero(vhdrWater);
vhdrWater.level[0] = 0;
vhdrWater.type[0] = Mdvx::VERT_TYPE_SURFACE;

MdvxField::setFieldName("WaterFlag", fhdrWater);
MdvxField::setFieldNameLong("Water, flag 1=water, 0=land, 2=missing", fhdrWater);
MdvxField::setUnits("", fhdrWater);
MdvxField *waterField = new MdvxField(fhdrWater, vhdrWater, water, false, false, false);
delete[] water;

// add water field to input object

_inMdvx.addField(waterField);

return 0;

Expand All @@ -429,7 +469,7 @@ int Ecco::_addTerrainHtField()
/////////////////////////////////////////////////////////
// add fields to the output object

void Ecco::_addFields()
void Ecco::_addFieldsToOutput()

{

Expand All @@ -445,8 +485,8 @@ void Ecco::_addFields()
fhdr2d.volume_size = planeSize32;
fhdr2d.encoding_type = Mdvx::ENCODING_FLOAT32;
fhdr2d.data_element_nbytes = 4;
fhdr2d.missing_data_value = _missing;
fhdr2d.bad_data_value = _missing;
fhdr2d.missing_data_value = _missingFl32;
fhdr2d.bad_data_value = _missingFl32;
fhdr2d.scale = 1.0;
fhdr2d.bias = 0.0;

Expand Down Expand Up @@ -555,8 +595,8 @@ void Ecco::_addFields()

Mdvx::field_header_t fhdr3d = dbzField->getFieldHeader();
Mdvx::vlevel_header_t vhdr3d = dbzField->getVlevelHeader();
fhdr3d.missing_data_value = _missing;
fhdr3d.bad_data_value = _missing;
fhdr3d.missing_data_value = _missingFl32;
fhdr3d.bad_data_value = _missingFl32;

if (_params.write_convective_dbz) {
// reflectivity only where convection has been identified
Expand Down Expand Up @@ -618,6 +658,26 @@ void Ecco::_addFields()
if (_params.write_clumping_debug_fields) {
_addClumpingDebugFields();
}

// add ht data if available

if (_params.use_terrain_ht_data) {
MdvxField *htFieldIn = _inMdvx.getField("TerrainHt");
if (htFieldIn != NULL) {
MdvxField *htFieldOut = new MdvxField(*htFieldIn);
htFieldOut->convertType(Mdvx::ENCODING_ASIS, Mdvx::COMPRESSION_GZIP);
_outMdvx.addField(htFieldOut);
}
// add water field
if (_params.add_water_layer) {
MdvxField *waterFieldIn = _inMdvx.getField("WaterFlag");
if (waterFieldIn != NULL) {
MdvxField *waterFieldOut = new MdvxField(*waterFieldIn);
waterFieldOut->convertType(Mdvx::ENCODING_ASIS, Mdvx::COMPRESSION_GZIP);
_outMdvx.addField(waterFieldOut);
}
}
}

}

Expand Down Expand Up @@ -708,8 +768,8 @@ void Ecco::_computeHts(double tempC,
fhdr.volume_size = planeSize32;
fhdr.encoding_type = Mdvx::ENCODING_FLOAT32;
fhdr.data_element_nbytes = 4;
fhdr.missing_data_value = _missing;
fhdr.bad_data_value = _missing;
fhdr.missing_data_value = _missingFl32;
fhdr.bad_data_value = _missingFl32;

Mdvx::vlevel_header_t vhdr;
MEM_zero(vhdr);
Expand Down Expand Up @@ -838,7 +898,7 @@ int Ecco::_doWrite()

// add fields

_addFields();
_addFieldsToOutput();

// write out

Expand Down Expand Up @@ -931,8 +991,8 @@ void Ecco::_addClumpingDebugFields()
fhdr2d.volume_size = planeSize32;
fhdr2d.encoding_type = Mdvx::ENCODING_FLOAT32;
fhdr2d.data_element_nbytes = 4;
fhdr2d.missing_data_value = _missing;
fhdr2d.bad_data_value = _missing;
fhdr2d.missing_data_value = _missingFl32;
fhdr2d.bad_data_value = _missingFl32;
fhdr2d.scale = 1.0;
fhdr2d.bias = 0.0;

Expand Down
4 changes: 2 additions & 2 deletions codebase/apps/radar/src/Ecco/Ecco.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ protected:

private:

static const fl32 _missing;
static const fl32 _missingFl32;

string _progName;
char *_paramsPath;
Expand All @@ -95,7 +95,7 @@ private:

int _doRead();
int _addTerrainHtField();
void _addFields();
void _addFieldsToOutput();
int _readTempProfile(time_t dbzTime,
const MdvxField *dbzField);
void _computeHts(double tempC,
Expand Down
36 changes: 24 additions & 12 deletions codebase/apps/radar/src/Ecco/Params.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -708,18 +708,6 @@
tt->single_val.s = tdrpStrDup("$(HOME)/data/terrain/DEM/uncompressed");
tt++;

// Parameter 'water_layer_dir'
// ctype is 'char*'

memset(tt, 0, sizeof(TDRPtable));
tt->ptype = STRING_TYPE;
tt->param_name = tdrpStrDup("water_layer_dir");
tt->descr = tdrpStrDup("Directory for water layer data in netCDF format.");
tt->help = tdrpStrDup("These files match the SRTM30 data set in spatial extent. The data resolution is 30 arc-seconds, or 120 per degree in lat/lon. The data is in bytes: 1 for water, 0 for not water.");
tt->val_offset = (char *) &water_layer_dir - &_start_;
tt->single_val.s = tdrpStrDup("$(HOME)/data/terrain/WATER");
tt++;

// Parameter 'check_adjacent_grid_cells'
// ctype is 'tdrp_bool_t'

Expand All @@ -744,6 +732,30 @@
tt->single_val.d = 1;
tt++;

// Parameter 'add_water_layer'
// ctype is 'tdrp_bool_t'

memset(tt, 0, sizeof(TDRPtable));
tt->ptype = BOOL_TYPE;
tt->param_name = tdrpStrDup("add_water_layer");
tt->descr = tdrpStrDup("Add water layer data to output file.");
tt->help = tdrpStrDup("If true, we read in the water layer for each grid location. This is a 1 or 0. It is 0 for land and 1 for water.");
tt->val_offset = (char *) &add_water_layer - &_start_;
tt->single_val.b = pTRUE;
tt++;

// Parameter 'water_layer_dir'
// ctype is 'char*'

memset(tt, 0, sizeof(TDRPtable));
tt->ptype = STRING_TYPE;
tt->param_name = tdrpStrDup("water_layer_dir");
tt->descr = tdrpStrDup("Directory for water layer data in netCDF format.");
tt->help = tdrpStrDup("These files match the SRTM30 data set in spatial extent. The data resolution is 30 arc-seconds, or 120 per degree in lat/lon. The data is in bytes: 1 for water, 0 for not water.");
tt->val_offset = (char *) &water_layer_dir - &_start_;
tt->single_val.s = tdrpStrDup("$(HOME)/data/terrain/WATER");
tt++;

// Parameter 'Comment 4'

memset(tt, 0, sizeof(TDRPtable));
Expand Down
8 changes: 5 additions & 3 deletions codebase/apps/radar/src/Ecco/Params.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -395,12 +395,14 @@ public:

char* srtm30_dem_dir;

char* water_layer_dir;

tdrp_bool_t check_adjacent_grid_cells;

double search_margin_km;

tdrp_bool_t add_water_layer;

char* water_layer_dir;

double min_valid_height;

double max_valid_height;
Expand Down Expand Up @@ -496,7 +498,7 @@ private:

void _init();

mutable TDRPtable _table[68];
mutable TDRPtable _table[69];

const char *_className;

Expand Down
8 changes: 6 additions & 2 deletions codebase/apps/radar/src/Ecco/SquareDegree.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -309,8 +309,12 @@ int SquareDegree::_readFromFile()

// read in water data

if (_readWaterFile(ulOffset)) {
return -1;
if (_params.add_water_layer) {
if (_readWaterFile(ulOffset)) {
return -1;
}
} else {
_waterAvail = false;
}

return 0;
Expand Down
Loading

0 comments on commit 000ee40

Please sign in to comment.