Skip to content

Commit

Permalink
Merge pull request #257 from KNMI/proj
Browse files Browse the repository at this point in the history 
Switch to new proj API
  • Loading branch information
@maartenplieger
maartenplieger committed May 11, 2023
2 parents b2331e7 + a41879e commit 755dcd9
Show file tree
Hide file tree
Showing 22 changed files with 274 additions and 616 deletions.
351 changes: 2 additions & 349 deletions CCDFDataModel/CCDFHDF5IO.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,8 @@

#include "CCDFHDF5IO.h"

#include <cmath>

const char *CDFHDF5Reader::className = "CDFHDF5Reader";

int CDFHDF5Reader::CustomForecastReader::readData(CDF::Variable *thisVar, size_t *start, size_t *count, ptrdiff_t *stride) {
Expand Down Expand Up @@ -365,8 +367,6 @@ int CDFHDF5Reader::open(const char *fileName) {
if (status == 1) return 1;
status = convertLSASAFtoCF();
if (status == 1) return 1;
status = convertKNMIH5VolScan();
if (status == 1) return 1;
}
if (b_EnableODIMHDF5toCFConversion) {
int status = convertODIMHDF5toCF();
Expand All @@ -381,10 +381,6 @@ int CDFHDF5Reader::close() {
delete forecastReader;
forecastReader = NULL;
}
if (volScanReader != NULL) {
delete volScanReader;
volScanReader = NULL;
}
fileIsOpen = false;
return 0;
}
Expand Down Expand Up @@ -1139,189 +1135,6 @@ int CDFHDF5Reader::convertLSASAFtoCF() {
return 0;
}

int CDFHDF5Reader::convertKNMIH5VolScan() {
return 0;
try {
/* Check if the image1.statistics variable and stat_cell_number is set */
CDF::Attribute *attr = cdfObject->getVariable("overview")->getAttribute("number_scan_groups");
int number_scan_groups;
attr->getData(&number_scan_groups, 1);
if (number_scan_groups == 0) return 0;
/* TODO: Return 1 in case it is not a point renderer*/
} catch (int e) {
return 0;
}

int scans[] = {8, 9, 10, 2, 11, 3, 12, 4, 13, 5, 14, 6, 7};
int nrscans = 13; // length(scans)
CT::string scan_params[] = {"KDP", "PhiDP", "RhoHV", "V", "W", "Z"};

for (size_t v = 0; v < cdfObject->variables.size(); v++) {
CDF::Variable *var = cdfObject->variables[v];
CT::string *terms = var->name.splitToArray(".");
if (terms->count > 1) {
if (terms[0].startsWith("scan") && terms[1].startsWith("scan_") && terms[1].endsWith("_data")) {
var->setAttributeText("ADAGUC_SKIP", "TRUE");
}
}
if (var->name.startsWith("visualisation")) {
var->setAttributeText("ADAGUC_SKIP", "TRUE");
}
}

size_t width = 700;
size_t height = 765;
double startX = 0;
double startY = -3649.980;
double step_x = 1.000;
double step_y = -1.000;

// For x
CDF::Dimension *dimX = new CDF::Dimension();
dimX->name = "x";
dimX->setSize(width);
cdfObject->addDimension(dimX);
CDF::Variable *varX = new CDF::Variable();
varX->setType(CDF_DOUBLE);
varX->name.copy("x");
varX->isDimension = true;
varX->dimensionlinks.push_back(dimX);
cdfObject->addVariable(varX);
CDF::allocateData(CDF_DOUBLE, &varX->data, dimX->length);
double *px = (double *)varX->data;
double xval = startX;
for (size_t i = 0; i < dimX->length; i++) {
*px++ = xval;
xval += step_x;
}
CDBDebug("X filled from %f to %f", startX, xval);

// For y
CDF::Dimension *dimY = new CDF::Dimension();
dimY->name = "y";
dimY->setSize(height);
cdfObject->addDimension(dimY);
CDF::Variable *varY = new CDF::Variable();
varY->setType(CDF_DOUBLE);
varY->name.copy("y");
varY->isDimension = true;
varY->dimensionlinks.push_back(dimY);
cdfObject->addVariable(varY);
CDF::allocateData(CDF_DOUBLE, &varY->data, dimY->length);
double *py = (double *)varY->data;
double yval = startY;
for (size_t i = 0; i < height; i++) {
*py++ = yval;
yval += step_y;
}
CDBDebug("Y filled from %f to %f", startY, yval);

// Create a new time dimension for the new 2D fields.
CDF::Dimension *dimT = new CDF::Dimension();
dimT->name = "time";
dimT->setSize(1);
cdfObject->addDimension(dimT);
CT::string time = cdfObject->getVariable("overview")->getAttribute("product_datetime_start")->getDataAsString();
CDF::Variable *timeVar = new CDF::Variable();
timeVar->setType(CDF_DOUBLE);
timeVar->name.copy("time");
timeVar->setAttributeText("standard_name", "time");
timeVar->setAttributeText("long_name", "time");
timeVar->isDimension = true;
char szStartTime[100];
CT::string time_units = "minutes since 2000-01-01 00:00:00\0";
CT::string h5Time = cdfObject->getVariable("overview")->getAttribute("product_datetime_start")->getDataAsString();
CDFHDF5Reader::HDF5ToADAGUCTime(szStartTime, h5Time.c_str());
// Set adaguc time
CTime ctime;
if (ctime.init(time_units, NULL) != 0) {
CDBError("Could not initialize CTIME: %s", time_units.c_str());
return 1;
}
double offset;
try {
offset = ctime.dateToOffset(ctime.stringToDate(szStartTime));
} catch (int e) {
CT::string message = CTime::getErrorMessage(e);
CDBError("CTime Exception %s", message.c_str());
return 1;
}
// CDBDebug("offset from %s = %f", szStartTime, offset);

timeVar->setAttributeText("units", time_units.c_str());
timeVar->dimensionlinks.push_back(dimT);
CDF::allocateData(CDF_DOUBLE, &timeVar->data, 1);
((double *)timeVar->data)[0] = offset;
cdfObject->addVariable(timeVar);

// Create a new elevation dimension for the new 2D fields.
CDF::Dimension *dimElevation = new CDF::Dimension();
dimElevation->name = "scan_elevation";
dimElevation->setSize(nrscans);
cdfObject->addDimension(dimElevation);
CDF::Variable *varElevation = new CDF::Variable();
varElevation->setType(CDF_DOUBLE);
varElevation->name.copy("scan_elevation");
CDF::Attribute *unit = new CDF::Attribute("units", "degrees");
varElevation->addAttribute(unit);
varElevation->isDimension = true;
varElevation->dimensionlinks.push_back(dimElevation);
cdfObject->addVariable(varElevation);
CDF::allocateData(CDF_DOUBLE, &varElevation->data, dimElevation->length);

CDF::Variable *varScan = new CDF::Variable();
varScan->setType(CDF_UINT);
varScan->name.copy("scan_number");
varScan->isDimension = false;
varScan->dimensionlinks.push_back(dimElevation);
cdfObject->addVariable(varScan);
CDF::allocateData(varScan->getType(), &varScan->data, dimElevation->length);

for (int i = 0; i < nrscans; i++) {
CT::string scanVarName;
scanVarName.print("scan%1d", scans[i]);
CDF::Variable *scanVar = cdfObject->getVariable(scanVarName.c_str());
float scanElevation;
scanVar->getAttribute("scan_elevation")->getData(&scanElevation, 1);
CDBDebug("scanelevation[%d] %d:%f", i, scans[i], scanElevation);
((double *)varElevation->data)[i] = scanElevation;
((unsigned int *)varScan->data)[i] = scans[i];
}

CDF::Variable *projection = new CDF::Variable();
projection->name = "projection";
projection->setType(CDF_CHAR);
projection->setAttributeText("projection_name", "STEREOGRAPHIC");
projection->setAttributeText("grid_mapping_name", "polar_stereographic");
projection->setAttributeText("proj4_params", KNMI_VOLSCAN_PROJ4);
cdfObject->addVariable(projection);

volScanReader = new CustomVolScanReader();
for (CT::string s : scan_params) {
CDBDebug("Adding variable %s", s.c_str());
CDF::Variable *var = new CDF::Variable();
var->setType(CDF_FLOAT);
;
var->name.copy(s);
cdfObject->addVariable(var);
var->setAttributeText("standard_name", s.c_str());
var->setAttributeText("long_name", s.c_str());
var->setAttributeText("grid_mapping", "projection");
var->setAttributeText("ADAGUC_VOL_SCAN", "TRUE");
var->setAttributeText("units", "fixedUnits");
float fillValue = FLT_MAX;
var->setAttribute("_FillValue", CDF_FLOAT, &fillValue, 1);
var->setCustomReader(volScanReader);

var->dimensionlinks.push_back(dimT);
var->dimensionlinks.push_back(dimElevation);
var->dimensionlinks.push_back(dimY);
var->dimensionlinks.push_back(dimX);
}

return 0;
}

int CDFHDF5Reader::readAttributes(std::vector<CDF::Attribute *> &attributes, hid_t HDF5_group) {
hid_t HDF5_attr_class, HDF5_attribute;
int dNumAttributes = H5Aget_num_attrs(HDF5_group);
Expand Down Expand Up @@ -1414,166 +1227,6 @@ int CDFHDF5Reader::readAttributes(std::vector<CDF::Attribute *> &attributes, hid
return 0;
}

int CDFHDF5Reader::CustomVolScanReader::getCalibrationParameters(CT::string formula, float &factor, float &offset) {
CDBDebug("Formula: %s", formula.c_str());
int rightPartFormulaPos = formula.indexOf("=");
int multiplicationSignPos = formula.indexOf("*");
int additionSignPos = formula.indexOf("+");
if (rightPartFormulaPos != -1 && multiplicationSignPos != -1 && additionSignPos != -1) {
factor = formula.substring(rightPartFormulaPos + 1, multiplicationSignPos).trim().toFloat();
offset = formula.substring(additionSignPos + 1, formula.length()).trim().toFloat();
return 0;
}
CDBDebug("Using default factor/offset from %s", formula.c_str());
factor = 1;
offset = 0;
return 1;
}

int CDFHDF5Reader::CustomVolScanReader::readData(CDF::Variable *thisVar, size_t *start, size_t *count, ptrdiff_t *stride) {
#ifdef CCDFHDF5IO_DEBUG
CDBDebug("READ data for %s called", thisVar->name.c_str());
#endif

size_t newstart[thisVar->dimensionlinks.size()];
for (size_t j = 0; j < thisVar->dimensionlinks.size(); j++) {
newstart[j] = start[j];
#ifdef CCDFHDF5IO_DEBUG
CDBDebug("%s %d %d %d %d", thisVar->dimensionlinks[j]->name.c_str(), j, start[j], count[j], stride[j]);
#endif
}
CDBDebug("Continuing for %s", thisVar->name.c_str());
newstart[0] = 0;
if (newstart[0] == 0) {
} // Silence TODO

CDF::Variable *varScanNumber = ((CDFObject *)thisVar->getParentCDFObject())->getVariable("scan_number");
unsigned int *scanNumbers = (unsigned int *)(varScanNumber->data);

CT::string varName;
varName.print("scan%d%sscan_%s_data", scanNumbers[(int)start[1]], CCDFHDF5IO_GROUPSEPARATOR, thisVar->name.c_str());
CDBDebug("varName: %s", varName.c_str());
CDF::Variable *var = ((CDFObject *)thisVar->getParentCDFObject())->getVariable(varName.c_str());

CDF::Dimension *xDim = ((CDFObject *)thisVar->getParentCDFObject())->getDimension("x");
size_t xSize = xDim->getSize();
CDF::Variable *xDimVar = ((CDFObject *)thisVar->getParentCDFObject())->getVariable("x");
double xMin = xDimVar->getDataAt<double>(0);
double xMax = xDimVar->getDataAt<double>(xSize - 1);
double xStep = (xMax - xMin) / (xSize - 1);

CDF::Dimension *yDim = ((CDFObject *)thisVar->getParentCDFObject())->getDimension("y");
size_t ySize = yDim->getSize();
CDF::Variable *yDimVar = ((CDFObject *)thisVar->getParentCDFObject())->getVariable("y");
double yMin = yDimVar->getDataAt<double>(0);
double yMax = yDimVar->getDataAt<double>(ySize - 1);
double yStep = (yMax - yMin) / (ySize - 1);

CT::string scanCalibrationVarName;
scanCalibrationVarName.print("scan%1d.calibration", scanNumbers[start[1]]);
CDF::Variable *scanCalibrationVar = ((CDFObject *)thisVar->getParentCDFObject())->getVariable(scanCalibrationVarName);
CT::string componentCalibrationStringName;
componentCalibrationStringName.print("calibration_%s_formulas", thisVar->name.c_str());
CT::string calibrationFormula = scanCalibrationVar->getAttribute(componentCalibrationStringName.c_str())->getDataAsString();
float factor, offset;
getCalibrationParameters(calibrationFormula, factor, offset);

CDBDebug("Using factor %f and offset %f for scan%d", factor, offset, scanNumbers[start[1]]);
float nodata = FLT_MAX;

CDBDebug("Start reading %s %d", var->name.c_str(), var->getType());

CDFType readType = var->getType();
CDBDebug("readType: %d", readType);
int status = var->readData(readType, false);
float *scanData = (float *)malloc(sizeof(float) * var->getSize());
float *pp = scanData;
for (size_t i = 0; i < var->getSize(); i++) {
if (((unsigned short *)var->data)[i] == 0) {
*pp++ = nodata;
} else {
*pp++ = ((unsigned short *)var->data)[i] * factor + offset;
}
}

CDBDebug("read and converted %d values (%d, %f)", var->getSize(), ((unsigned short *)(var->data))[0], *scanData);
// int status = var->readData(readType, newstart, count, stride, true);

if (status != 0) {
CDBError("CustomVolScanReader: Unable to read variable %s", thisVar->name.c_str());
return 1;
}

CDF::freeData(&thisVar->data);
int size = 1;
for (size_t j = 0; j < thisVar->dimensionlinks.size(); j++) {
size *= int((float(count[j]) / float(stride[j])) + 0.5);
}
thisVar->setSize(size);

CDBDebug("Reading size: %d", size);

CDF::allocateData(thisVar->getType(), &thisVar->data, size);

float radarLonLat[2];
((CDFObject *)thisVar->getParentCDFObject())->getVariable("radar1")->getAttribute("radar_location")->getData<float>(radarLonLat, 2);
float radarLon = radarLonLat[0];
float radarLat = radarLonLat[1];

CT::string scanName;
scanName.print("scan%1d", scanNumbers[start[1]]);
CDF::Variable *scanVar = ((CDFObject *)thisVar->getParentCDFObject())->getVariable(scanName);
float scan_elevation;
scanVar->getAttribute("scan_elevation")->getData<float>(&scan_elevation, 1);
int scan_nrang;
scanVar->getAttribute("scan_number_range")->getData<int>(&scan_nrang, 1);
int scan_nazim;
scanVar->getAttribute("scan_number_azim")->getData<int>(&scan_nazim, 1);
float scan_rscale;
scanVar->getAttribute("scan_range_bin")->getData<float>(&scan_rscale, 1);
float scan_ascale;
scanVar->getAttribute("scan_azim_bin")->getData<float>(&scan_ascale, 1);
CDBDebug("scanparameters [%d:%d] %f, %d, %d, %f, %f", start[1], scanNumbers[start[1]], scan_elevation, scan_nrang, scan_nazim, scan_rscale, scan_ascale);

/*Setting geographical projection parameters of input Cartesian grid.*/
CT::string scanProj4;
scanProj4.print("+proj=aeqd +a=6378.137 +b=6356.752 +R_A +lat_0=%.3f +lon_0=%.3f +x_0=0 +y_0=0", radarLat, radarLon);
projPJ projAeqd = pj_init_plus(scanProj4.c_str());
if (projAeqd == NULL) {
}

CT::string compositeProj4 = KNMI_VOLSCAN_PROJ4;
projPJ projComposite = pj_init_plus(compositeProj4.c_str());

double x, y;
float rang, azim;
int ir, ia;
float *p = (float *)thisVar->data;
float *pScan = scanData;

CDBDebug("Filling %d * %d cells", count[2] * stride[2], count[3] * stride[3]);
for (size_t row = start[2]; row < count[2]; row += stride[2]) {
for (size_t col = start[3]; col < count[3]; col += stride[3]) {
x = xMin + (col + 0.5) * xStep;
y = yMin + (row + 0.5) * yStep;
pj_transform(projComposite, projAeqd, 1, 1, &x, &y, NULL);
rang = sqrt(x * x + y * y);
azim = atan2(x, y) * 180.0 / M_PI;
ir = (int)(rang / scan_rscale);
ia = (int)(azim / scan_ascale);
ia = (ia + scan_nazim) % scan_nazim;
if (ir < scan_nrang)
*p++ = pScan[ir + ia * scan_nrang];
else
*p++ = nodata;
}
}

free(scanData);

return 0;
}

int CDFHDF5Reader::convertKNMIHDF5toCF() {

CDF::Variable *geo = cdfObject->getVariableNE("geographic");
Expand Down
Loading

0 comments on commit 755dcd9

Please sign in to comment.