diff --git a/Dockerfile b/Dockerfile index 34f0e6e51..7c5615df4 100755 --- a/Dockerfile +++ b/Dockerfile @@ -6,7 +6,7 @@ USER root LABEL maintainer="adaguc@knmi.nl" # Version should be same as in Definitions.h -LABEL version="2.28.3" +LABEL version="2.28.4" # Try to update image packages RUN apt-get -q -y update \ diff --git a/NEWS.md b/NEWS.md index 6fd96ce54..3f63095b1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,6 @@ +**Version 2.28.4 2024-11-04** +- Added Solar Terminator post-processor to use in combination with LiveUpdate layer type. + **Version 2.28.3 2024-10-24** - Fix bug where directory reader could not figure out the file type (opendir ent->d_type == DT_UNKNOWN) diff --git a/adagucserverEC/CDataPostProcessors/CDataPostProcessor.cpp b/adagucserverEC/CDataPostProcessors/CDataPostProcessor.cpp index 26e4ccaf3..d328ad482 100644 --- a/adagucserverEC/CDataPostProcessors/CDataPostProcessor.cpp +++ b/adagucserverEC/CDataPostProcessors/CDataPostProcessor.cpp @@ -12,6 +12,7 @@ #include "CDataPostProcessors_MSGCPP.h" #include "CDataPostProcessor_CDPDBZtoRR.h" #include "CDataPostProcessor_AddFeatures.h" +#include "CDataPostProcessor_SolarTerminator.h" extern CDPPExecutor cdppExecutorInstance; CDPPExecutor cdppExecutorInstance; @@ -37,6 +38,7 @@ CDPPExecutor::CDPPExecutor() { dataPostProcessorList->push_back(new CDPPOperator()); dataPostProcessorList->push_back(new CDPPWFP()); dataPostProcessorList->push_back(new CDPPWindSpeedKnotsToMs()); + dataPostProcessorList->push_back(new CDPPSolarTerminator()); } CDPPExecutor::~CDPPExecutor() { diff --git a/adagucserverEC/CDataPostProcessors/CDataPostProcessor_SolarTerminator.cpp b/adagucserverEC/CDataPostProcessors/CDataPostProcessor_SolarTerminator.cpp new file mode 100644 index 000000000..c2ae55b03 --- /dev/null +++ b/adagucserverEC/CDataPostProcessors/CDataPostProcessor_SolarTerminator.cpp @@ -0,0 +1,219 @@ +#include "CDataPostProcessor_SolarTerminator.h" +#include "solar/solar_terminator.h" +#include +#include +#include "CTime.h" +#include "CImageWarper.h" + +#include +#include +#include +#include +#include + +#include + +/************************/ +/* CDPPSolarTerminator */ +/************************/ +const char *CDPPSolarTerminator::className = "CDPPSolarTerminator"; + +const char *CDPPSolarTerminator::getId() { return "solarterminator"; } + +int CDPPSolarTerminator::isApplicable(CServerConfig::XMLE_DataPostProc *proc, CDataSource *dataSource, int mode) { + if (proc->attr.algorithm.equals("solarterminator")) { + if (dataSource->getNumDataObjects() < 1 && mode == CDATAPOSTPROCESSOR_RUNBEFOREREADING) { + CDBError("1 variable is needed for solarterminator, found %d", dataSource->getNumDataObjects()); + return CDATAPOSTPROCESSOR_CONSTRAINTSNOTMET; + } + return CDATAPOSTPROCESSOR_RUNAFTERREADING | CDATAPOSTPROCESSOR_RUNBEFOREREADING; + } + return CDATAPOSTPROCESSOR_NOTAPPLICABLE; +} + +int CDPPSolarTerminator::execute(CServerConfig::XMLE_DataPostProc *proc, CDataSource *dataSource, int mode) { + if ((isApplicable(proc, dataSource, mode) & mode) == false) { + return -1; + } + + double currentOffset = 1; + + if (!dataSource->srvParams->requestDims.empty()) { + CT::string timestampStr = dataSource->srvParams->requestDims[0]->value.c_str(); + currentOffset = CTime::getEpochTimeFromDateString(dataSource->srvParams->requestDims[0]->value); + } + + if (mode == CDATAPOSTPROCESSOR_RUNBEFOREREADING) { + CT::string newVariableName = "SolT"; + + CDataSource::DataObject *newDataObject = dataSource->getDataObject(0); + newDataObject->variableName.copy(newVariableName.c_str()); + + // Copy bounding box of screen + auto *geo = dataSource->srvParams->Geo; + double dfBBOX[] = {geo->dfBBOX[0], geo->dfBBOX[1], geo->dfBBOX[2], geo->dfBBOX[3]}; + size_t width = geo->dWidth; + size_t height = geo->dHeight; + + dataSource->nativeProj4 = geo->CRS; + dataSource->dWidth = geo->dWidth; + dataSource->dHeight = geo->dHeight; + dataSource->dfBBOX[0] = geo->dfBBOX[0]; + dataSource->dfBBOX[1] = geo->dfBBOX[1]; + dataSource->dfBBOX[2] = geo->dfBBOX[2]; + dataSource->dfBBOX[3] = geo->dfBBOX[3]; + + // Create new dimensions and variables (X,Y,T) + + CDF::Dimension *dimX = new CDF::Dimension(); + dimX->name = "xet"; + dimX->setSize(width); + newDataObject->cdfObject->addDimension(dimX); + + // Define the X variable using the X dimension + CDF::Variable *varX = new CDF::Variable(); + varX->setType(CDF_DOUBLE); + varX->name.copy("xet"); + varX->isDimension = true; + varX->dimensionlinks.push_back(dimX); + newDataObject->cdfObject->addVariable(varX); + CDF::allocateData(CDF_DOUBLE, &varX->data, dimX->length); + + // Set the bbox in the data, since the virtual grid is 2x2 pixels we can directly apply the bbox + ((double *)varX->data)[0] = dfBBOX[0]; + ((double *)varX->data)[1] = dfBBOX[2]; + + // For y dimension + CDF::Dimension *dimY = new CDF::Dimension(); + dimY->name = "yet"; + dimY->setSize(height); + newDataObject->cdfObject->addDimension(dimY); + + // Define the Y variable using the X dimension + CDF::Variable *varY = new CDF::Variable(); + varY->setType(CDF_DOUBLE); + varY->name.copy("yet"); + varY->isDimension = true; + varY->dimensionlinks.push_back(dimY); + newDataObject->cdfObject->addVariable(varY); + CDF::allocateData(CDF_DOUBLE, &varY->data, dimY->length); + + ((double *)varY->data)[0] = dfBBOX[1]; + ((double *)varY->data)[1] = dfBBOX[3]; + + // For time dimension + CDF::Dimension *dimTime = new CDF::Dimension(); + dimTime->name = "time"; + dimTime->setSize(10); // 24 * 6); // Last day every 10 minutes + newDataObject->cdfObject->addDimension(dimTime); + + // Define the Y variable using the X dimension + CDF::Variable *varTime = new CDF::Variable(); + varTime->setType(CDF_DOUBLE); + varTime->name.copy("time"); + varTime->isDimension = true; + varTime->dimensionlinks.push_back(dimTime); + newDataObject->cdfObject->addVariable(varTime); + varTime->setAttributeText("units", "seconds since 1970"); + CTime epochCTime; + epochCTime.init("seconds since 1970-01-01 0:0:0", NULL); + CDF::allocateData(CDF_DOUBLE, &varTime->data, dimTime->length); + + for (int off = 0; off < 10; off++) { + // Every 10 minutes for a day + double timestep = epochCTime.quantizeTimeToISO8601(currentOffset - off * 60 * 10, "PT30M", "low"); + ((double *)varTime->data)[off] = timestep; // timestep; + } + + dataSource->getDataObject(0)->cdfVariable->dimensionlinks.push_back(dimTime); + // Define the Solar Terminator variable using the defined dimensions, and set the right attributes + CDF::Variable *solTVar = new CDF::Variable(); + solTVar->setType(CDF_FLOAT); + float fillValue[] = {-1}; + solTVar->setAttribute("_FillValue", solTVar->getType(), fillValue, 1); + solTVar->dimensionlinks.push_back(dimY); + solTVar->dimensionlinks.push_back(dimX); + solTVar->setType(CDF_FLOAT); + solTVar->name = "SolT"; + CDBDebug("Setting units"); + solTVar->setAttributeText("units", "categories"); + solTVar->setAttributeText("grid_mapping", "projection"); + newDataObject->cdfObject->addVariable(solTVar); + newDataObject->cdfVariable = solTVar; + + newDataObject->cdfVariable->setCustomReader(CDF::Variable::CustomMemoryReaderInstance); + + newDataObject->cdfVariable->setSize(dataSource->dWidth * dataSource->dHeight); + + // Make the width and height of the new 2D adaguc field the same as the viewing window + dataSource->dWidth = dataSource->srvParams->Geo->dWidth; + dataSource->dHeight = dataSource->srvParams->Geo->dHeight; + + // Width and height of the dataSource need to be at least 2 in this case. + if (dataSource->dWidth < 2) dataSource->dWidth = 2; + if (dataSource->dHeight < 2) dataSource->dHeight = 2; + + // Get the X and Y dimensions previousely defined and adjust them to the new settings and new grid (Grid in screenview space) + dimX->setSize(dataSource->dWidth); + dimY->setSize(dataSource->dHeight); + + // Re-allocate data for these coordinate variables with the new grid size + CDF::allocateData(CDF_DOUBLE, &varX->data, dimX->length); + CDF::allocateData(CDF_DOUBLE, &varY->data, dimY->length); + + // Calculate the gridsize, allocate data and fill the data with a fillvalue + size_t fieldSize = dimX->length * dimY->length; + newDataObject->cdfVariable->setSize(fieldSize); + CDF::allocateData(newDataObject->cdfVariable->getType(), &(newDataObject->cdfVariable->data), fieldSize); + CDF::fill(newDataObject->cdfVariable->data, newDataObject->cdfVariable->getType(), fillValue[0], fieldSize); + + // Calculate cellsize and offset of the echo toppen (ET) 2D virtual grid, using the same grid as the screenspace + double cellSizeX = (dataSource->srvParams->Geo->dfBBOX[2] - dataSource->srvParams->Geo->dfBBOX[0]) / double(dataSource->dWidth); + double cellSizeY = (dataSource->srvParams->Geo->dfBBOX[3] - dataSource->srvParams->Geo->dfBBOX[1]) / double(dataSource->dHeight); + double offsetX = dataSource->srvParams->Geo->dfBBOX[0]; + double offsetY = dataSource->srvParams->Geo->dfBBOX[1]; + + // Fill in the X and Y dimensions with the array of coordinates + for (size_t j = 0; j < dimX->length; j++) { + double x = offsetX + double(j) * cellSizeX + cellSizeX / 2; + ((double *)varX->data)[j] = x; + } + for (size_t j = 0; j < dimY->length; j++) { + double y = offsetY + double(j) * cellSizeY + cellSizeY / 2; + ((double *)varY->data)[j] = y; + } + } + if (mode == CDATAPOSTPROCESSOR_RUNAFTERREADING) { + CDBDebug("CDATAPOSTPROCESSOR_RUNAFTERREADING::Applying SOLARTERMINATOR"); + size_t l = (size_t)dataSource->dHeight * (size_t)dataSource->dWidth; + CDF::allocateData(dataSource->getDataObject(0)->cdfVariable->getType(), &dataSource->getDataObject(0)->cdfVariable->data, l); + + float *result = (float *)dataSource->getDataObject(0)->cdfVariable->data; + + CImageWarper imageWarper; + int status = imageWarper.initreproj(dataSource, dataSource->srvParams->Geo, &dataSource->srvParams->cfg->Projection); + if (status != 0) { + CDBError("Unable to init projection"); + return 1; + } + + for (size_t j = 0; j < l; j++) { + int px = j % dataSource->dWidth; + int py = j / dataSource->dWidth; + + double lonRange = dataSource->dfBBOX[2] - dataSource->dfBBOX[0]; + double latRange = dataSource->dfBBOX[1] - dataSource->dfBBOX[3]; + + // Projection coordinates (works in EPSG 4326) + double geox = (lonRange / dataSource->dWidth) * px + dataSource->dfBBOX[0]; + double geoy = (latRange / dataSource->dHeight) * py + dataSource->dfBBOX[3]; + + // Transform EPG:3857 coordinates into latlon + imageWarper.reprojToLatLon(geox, geoy); + + // Select final value based on solar zenith angle + result[j] = static_cast(getDayTimeCategory(getSolarZenithAngle(geoy, geox, currentOffset))); + } + } + return 0; +} \ No newline at end of file diff --git a/adagucserverEC/CDataPostProcessors/CDataPostProcessor_SolarTerminator.h b/adagucserverEC/CDataPostProcessors/CDataPostProcessor_SolarTerminator.h new file mode 100644 index 000000000..937a85b23 --- /dev/null +++ b/adagucserverEC/CDataPostProcessors/CDataPostProcessor_SolarTerminator.h @@ -0,0 +1,20 @@ +#include "CDataPostProcessor.h" + +#ifndef CDATAPOSTPROCESSOR_SOLARTERMINATOR_H +#define CDATAPOSTPROCESSOR_SOLARTERMINATOR_H +/** + * SolarTerminator algorithm + */ + +class CDPPSolarTerminator : public CDPPInterface { +private: + DEF_ERRORFUNCTION(); + +public: + virtual const char *getId(); + virtual int isApplicable(CServerConfig::XMLE_DataPostProc *proc, CDataSource *dataSource, int mode); + virtual int execute(CServerConfig::XMLE_DataPostProc *proc, CDataSource *dataSource, int mode); + virtual int execute(CServerConfig::XMLE_DataPostProc *, CDataSource *, int, double *, size_t) { return 1; } // TODO: Still need to implement +}; + +#endif \ No newline at end of file diff --git a/adagucserverEC/CDataPostProcessors/solar/solar_terminator.cpp b/adagucserverEC/CDataPostProcessors/solar/solar_terminator.cpp new file mode 100644 index 000000000..cb9432ff8 --- /dev/null +++ b/adagucserverEC/CDataPostProcessors/solar/solar_terminator.cpp @@ -0,0 +1,167 @@ +#include +#include +#include +#include + +int getDayTimeCategory(double zenithAngle) { + if (zenithAngle < 90) + return 4; // "Daylight"; + else if (zenithAngle < 96) + return 3; // "Civil twilight"; + else if (zenithAngle < 102) + return 2; // "Nautical twilight"; + else if (zenithAngle < 108) + return 1; // "Astronomical twilight"; + else + return 0; // "Night"; +} + +const double DEG_TO_RAD = M_PI / 180.0; +const double RAD_TO_DEG = 180.0 / M_PI; + +double get_julian_day(double timestamp) { + // Transform from UNIX timestamp to Julian Day + // January 1, 1970 - Julian day no. 2440587.5 + double jd = 2440587.5 + timestamp / 86400.0; + return jd; +} + +double get_julian_century(double jd) { + // One Julian Century = 36525 days + // January 1, 2000 = Julian day no. 2451545 + // Accuracy here is of utmost importance (T is expressed in centuries) + double T = (jd - 2451545.0) / 36525.0; + return T; +} + +double get_solar_mean_longitude(double T) { + // Measured in degrees + double L0 = std::fmod(280.46646 + T * (36000.76983 + 0.0003032 * T), 360); + while (L0 > 360.0) { + L0 -= 360.0; + } + while (L0 < 0.0) { + L0 += 360.0; + } + return L0; +} + +double get_solar_mean_anomaly(double T) { + // Measured in degrees + double M = 357.52910 + 35999.05030 * T - 0.0001559 * T * T - 0.00000048 * T * T * T; + return M; +} + +double get_solar_equation_of_center(double T, double M) { + // Measured in degrees + double M_rad = M * DEG_TO_RAD; + double C = (1.914602 - 0.004817 * T - 0.000014 * T * T) * sin(M_rad) + (0.019993 - 0.000101 * T) * sin(2.0 * M_rad) + 0.000290 * sin(3.0 * M_rad); + return C; +} + +double get_long_asc_node(double T) { + // Measured in degrees + double omega = 125.04 - 1934.136 * T; + return omega; +} + +double get_solar_app_longitude(double L, double omega) { + // Measured in degrees + double Lapp = L - 0.00569 - 0.00478 * sin(omega * DEG_TO_RAD); + return Lapp; +} + +double get_solar_mean_obliquity_ecliptic(double T) { + // Measured in degrees + double e0 = 23 + 26 / 60.0 + 21.448 / 3600 - 46.8150 * T / 3600 - 0.00059 * T * T / 3600 + 0.00183 * T * T * T / 3600 + 0.00256 * cos((125.04 - 1934.136 * T) * DEG_TO_RAD); + return e0; +} + +double get_solar_obliquity_correction(double e0, double omega) { + // Measured in degrees + double e = e0 + 0.00256 * cos(omega * DEG_TO_RAD); + return e; +} + +double get_solar_declination(double e, double Lapp) { + // Measured in degrees + double delta = (asin(sin(e * DEG_TO_RAD) * sin(Lapp * DEG_TO_RAD))) * RAD_TO_DEG; + return delta; +} + +double get_right_ascension(double L, double e) { + double ra = fmod(atan2(sin(L * DEG_TO_RAD) * cos(e * DEG_TO_RAD), cos(L * DEG_TO_RAD)) * RAD_TO_DEG, 360.0); + return ra; +} + +double get_greenwich_hour_angle(double jd, double T) { + double gha = 280.46061837 + 360.98564736629 * (jd - 2451545) + 0.000387933 * T * T - T * T * T / 38710000; + return gha; +} + +double get_local_hour_angle(double gha, double longitude, double ra) { + double lha = gha + longitude - ra; + lha = fmod(lha, 360.0); + if (lha < 0.0) { + lha += 360.0; // Ensure positive LHA + } + return lha; +} + +double getSolarZenithAngle(double lat, double lon, double timestamp) { + // Using formula cos(θ) = sin(δ) × sin(φ) + cos(δ) × cos(φ) × cos(H) + // θ is the solar zenith angle, + // δ is the solar declination + // φ is the observer's latitude, + // H is the hour angle of the Sun. + // Get current Julian Day + double jd = get_julian_day(timestamp); + + // Get current Julian Century + double T = get_julian_century(jd); + + // Get Mean Longitude of the Sun + double L0 = get_solar_mean_longitude(T); + + // Get Mean Anomaly of the Sun (page 151) + double M = get_solar_mean_anomaly(T); + + // Get Equation of Center of the Sun + double C = get_solar_equation_of_center(T, M); + + // Get True Longitude of the Sun (λ) (page 152) + double L = L0 + C; + + // Get Longitude of the Ascending Node + double omega = get_long_asc_node(T); + + // Get Apparent Longitude of the Sun + double Lapp = get_solar_app_longitude(L, omega); + + // Get Mean Obliquity of the Ecliptic of the Sun + double e0 = get_solar_mean_obliquity_ecliptic(T); + + // Get Obliquity correction + double e = get_solar_obliquity_correction(e0, omega); + + // Get Declination of the Sun + double delta = get_solar_declination(e, Lapp); + + // Get Greenwich Hour Angle + double gha = get_greenwich_hour_angle(jd, T); + + // Get Right Ascension Angle + double ra = get_right_ascension(L, e); + + // Get Local Hour Angle + double lha = get_local_hour_angle(gha, lon, ra); + + // Final calculation of solar zenith angle + // cos(θ) = sin(δ) × sin(φ) + cos(δ) × cos(φ) × cos(H) + // θ is the solar zenith angle - called theta here + // δ is the solar declination - called delta here + // φ is the observer's latitude - called lat here + // H is the hour angle of the Sun - called lha here + double theta = RAD_TO_DEG * acos(sin(delta * DEG_TO_RAD) * sin(lat * DEG_TO_RAD) + cos(delta * DEG_TO_RAD) * cos(lat * DEG_TO_RAD) * cos(lha * DEG_TO_RAD)); + return theta; +} \ No newline at end of file diff --git a/adagucserverEC/CDataPostProcessors/solar/solar_terminator.h b/adagucserverEC/CDataPostProcessors/solar/solar_terminator.h new file mode 100644 index 000000000..b925c8b24 --- /dev/null +++ b/adagucserverEC/CDataPostProcessors/solar/solar_terminator.h @@ -0,0 +1,7 @@ +#ifndef SOLAR_POSITION_H +#define SOLAR_POSITION_H + +double getSolarZenithAngle(double lat, double lon, double timestamp); +int getDayTimeCategory(double zenithAngle); + +#endif // SOLAR_POSITION_H diff --git a/adagucserverEC/CMakeLists.txt b/adagucserverEC/CMakeLists.txt index 69b234c70..a5bb10147 100644 --- a/adagucserverEC/CMakeLists.txt +++ b/adagucserverEC/CMakeLists.txt @@ -106,6 +106,10 @@ add_library( CDataPostProcessors/CDataPostProcessor_IncludeLayer.h CDataPostProcessors/CDataPostProcessor_Operator.cpp CDataPostProcessors/CDataPostProcessor_Operator.h + CDataPostProcessors/CDataPostProcessor_SolarTerminator.h + CDataPostProcessors/CDataPostProcessor_SolarTerminator.cpp + CDataPostProcessors/solar/solar_terminator.h + CDataPostProcessors/solar/solar_terminator.cpp CDataPostProcessors/CDataPostProcessor_WFP.h CDataPostProcessors/CDataPostProcessor_WFP.cpp CDataReader.cpp diff --git a/adagucserverEC/CRequest.cpp b/adagucserverEC/CRequest.cpp index 8c3d40fe1..0bdb33a53 100644 --- a/adagucserverEC/CRequest.cpp +++ b/adagucserverEC/CRequest.cpp @@ -1676,6 +1676,7 @@ int CRequest::process_all_layers() { // dataSources[j]->getCDFDims()->addDimension("none","0",0); } if (dataSources[j]->dLayerType == CConfigReaderLayerTypeLiveUpdate) { + // This layer has no dimensions, but we need to add one timestep with data in order to make the next code work. layerTypeLiveUpdateConfigureDimensionsInDataSource(dataSources[j]); } } @@ -1715,7 +1716,7 @@ int CRequest::process_all_layers() { /* Handle WMS Getmap database request */ /**************************************/ if (dataSources[j]->dLayerType == CConfigReaderLayerTypeDataBase || dataSources[j]->dLayerType == CConfigReaderLayerTypeStyled || dataSources[j]->dLayerType == CConfigReaderLayerTypeCascaded || - dataSources[j]->dLayerType == CConfigReaderLayerTypeBaseLayer) { + dataSources[j]->dLayerType == CConfigReaderLayerTypeBaseLayer || (dataSources[j]->dLayerType == CConfigReaderLayerTypeLiveUpdate && srvParam->requestType != REQUEST_WMS_GETMAP)) { try { for (size_t d = 0; d < dataSources.size(); d++) { dataSources[d]->setTimeStep(0); @@ -2161,11 +2162,7 @@ int CRequest::process_all_layers() { return 1; } } else if (dataSources[j]->dLayerType == CConfigReaderLayerTypeLiveUpdate) { - // Render the current time in an image for testing purpose / frontend development - CDrawImage image; - layerTypeLiveUpdateRenderIntoDrawImage(&image, srvParam); - printf("%s%c%c\n", "Content-Type:image/png", 13, 10); - image.printImagePng8(true); + layerTypeLiveUpdateRender(dataSources[j], srvParam); } else { CDBError("Unknown layer type"); } diff --git a/adagucserverEC/CXMLGen.cpp b/adagucserverEC/CXMLGen.cpp index ebd21a808..c6f2752aa 100644 --- a/adagucserverEC/CXMLGen.cpp +++ b/adagucserverEC/CXMLGen.cpp @@ -215,7 +215,8 @@ int CXMLGen::getWMS_1_1_1_Capabilities(CT::string *XMLDoc, std::vectorhasError == 0) { - XMLDoc->printconcat("\n", layer->layerMetadata.isQueryable, layer->dataSource->dLayerType == CConfigReaderLayerTypeCascaded ? 1 : 0); + XMLDoc->printconcat("\n", layer->layerMetadata.isQueryable, + layer->dataSource->dLayerType == CConfigReaderLayerTypeCascaded && layer->dataSource->dLayerType == CConfigReaderLayerTypeLiveUpdate ? 1 : 0); XMLDoc->concat(""); XMLDoc->concat(&layer->layerMetadata.name); XMLDoc->concat("\n"); @@ -585,7 +586,8 @@ int CXMLGen::getWMS_1_3_0_Capabilities(CT::string *XMLDoc, std::vectorhasError == 0) { - XMLDoc->printconcat("\n", layer->layerMetadata.isQueryable, layer->dataSource->dLayerType == CConfigReaderLayerTypeCascaded ? 1 : 0); + XMLDoc->printconcat("\n", layer->layerMetadata.isQueryable, + layer->dataSource->dLayerType == CConfigReaderLayerTypeCascaded && layer->dataSource->dLayerType == CConfigReaderLayerTypeLiveUpdate ? 1 : 0); XMLDoc->concat(""); XMLDoc->concat(&layer->layerMetadata.name); XMLDoc->concat("\n"); diff --git a/adagucserverEC/Definitions.h b/adagucserverEC/Definitions.h index 2b37dddac..5b3369b3b 100755 --- a/adagucserverEC/Definitions.h +++ b/adagucserverEC/Definitions.h @@ -28,7 +28,7 @@ #ifndef Definitions_H #define Definitions_H -#define ADAGUCSERVER_VERSION "2.28.3" // Please also update in the Dockerfile to the same version +#define ADAGUCSERVER_VERSION "2.28.4" // Please also update in the Dockerfile to the same version // CConfigReaderLayerType #define CConfigReaderLayerTypeUnknown 0 diff --git a/adagucserverEC/LayerTypeLiveUpdate/LayerTypeLiveUpdate.cpp b/adagucserverEC/LayerTypeLiveUpdate/LayerTypeLiveUpdate.cpp index 9143f3342..a565906ff 100644 --- a/adagucserverEC/LayerTypeLiveUpdate/LayerTypeLiveUpdate.cpp +++ b/adagucserverEC/LayerTypeLiveUpdate/LayerTypeLiveUpdate.cpp @@ -1,4 +1,5 @@ #include "LayerTypeLiveUpdate.h" +#include "CServerParams.h" int layerTypeLiveUpdateConfigureDimensionsInDataSource(CDataSource *dataSource) { // This layer has no dimensions, but we need to add one timestep with data in order to make the next code work. @@ -14,11 +15,49 @@ int layerTypeLiveUpdateConfigureDimensionsInDataSource(CDataSource *dataSource) requiredDim->value = "2020-01-02T00:00:00Z"; dataSource->requiredDims.push_back(requiredDim); } - dataSource->addStep("", NULL); - dataSource->getCDFDims()->addDimension("none", "0", 0); + // Basic case of liveupdate layer + if (dataSource->cfgLayer->DataPostProc.empty()) { + // // Add step to empty file + dataSource->addStep("", NULL); + dataSource->getCDFDims()->addDimension("none", "0", 0); + } else { + // Case of liveupdate layers with data post processors (such as solar terminator) + // Currently one (dummy) file is required + std::vector fileList; + if (!dataSource->cfgLayer->FilePath.empty()) { + try { + fileList = CDBFileScanner::searchFileNames(dataSource->cfgLayer->FilePath[0]->value.c_str(), dataSource->cfgLayer->FilePath[0]->attr.filter, NULL); + } catch (int e) { + CDBError("Could not find any filename"); + return 1; + } + + if (fileList.size() == 0) { + CDBError("fileList.size()==0"); + return 1; + } + dataSource->addStep(fileList[0].c_str(), NULL); + } + } return 0; } +int layerTypeLiveUpdateRender(CDataSource *dataSource, CServerParams *srvParam) { + CDBDebug("in special liveupdate case with timesteps %d", dataSource->getNumTimeSteps()); + + if (dataSource->cfgLayer->DataPostProc.empty()) { + // Demo case: render the current time in an image for testing purposes / frontend development + CDrawImage image; + layerTypeLiveUpdateRenderIntoDrawImage(&image, srvParam); + printf("%s%c%c\n", "Content-Type:image/png", 13, 10); + CDBDebug("***Number of timesteps %d", dataSource->getNumTimeSteps()); + return image.printImagePng8(true); + } else { + // General case: Liveupdate with some data postprocessors + return layerTypeLiveUpdateRenderIntoImageDataWriter(dataSource, srvParam); + } +} + int layerTypeLiveUpdateRenderIntoDrawImage(CDrawImage *image, CServerParams *srvParam) { image->enableTransparency(true); image->setTrueColor(true); @@ -39,6 +78,26 @@ int layerTypeLiveUpdateRenderIntoDrawImage(CDrawImage *image, CServerParams *srv return 0; } +int layerTypeLiveUpdateRenderIntoImageDataWriter(CDataSource *dataSource, CServerParams *srvParam) { + // General case: Liveupdate with some data postprocessors + // Covers case of Solar Terminator + CImageDataWriter imageDataWriter; + int status = imageDataWriter.init(srvParam, dataSource, dataSource->getNumTimeSteps()); + CDBDebug("Init imageDataWriter status %d", status); + + if (dataSource->getNumTimeSteps() > 1) { + CDBDebug("Status from create animation was %d", imageDataWriter.createAnimation()); + } + + std::vector dataSourceRef = {dataSource}; + status = imageDataWriter.addData(dataSourceRef); + CDBDebug("Adding data status was %d", status); + + status = imageDataWriter.end(); + CDBDebug("Ending image data writing with status %d", status); + return status; +} + int layerTypeLiveUpdateConfigureWMSLayerForGetCapabilities(MetadataLayer *metadataLayer) { if (metadataLayer->dataSource->cfgLayer->Title.size() != 0) { metadataLayer->layerMetadata.title.copy(metadataLayer->dataSource->cfgLayer->Title[0]->value.c_str()); @@ -49,11 +108,11 @@ int layerTypeLiveUpdateConfigureWMSLayerForGetCapabilities(MetadataLayer *metada timeInstance.init("seconds since 1970", "standard"); double epochTime = timeInstance.getEpochTimeFromDateString(CTime::currentDateTime()); // CTime::Date cdate = timeInstance.getDate(epochTime); - double startTimeOffset = timeInstance.quantizeTimeToISO8601(epochTime - 3600, "PT1S", "low"); - double stopTimeOffset = timeInstance.quantizeTimeToISO8601(epochTime, "PT1S", "low"); + double startTimeOffset = timeInstance.quantizeTimeToISO8601(epochTime - 3600 * 24 * 365, "PT10M", "low"); + double stopTimeOffset = timeInstance.quantizeTimeToISO8601(epochTime, "PT10M", "low"); CT::string startTime = timeInstance.dateToISOString(timeInstance.offsetToDate(startTimeOffset)); CT::string stopTime = timeInstance.dateToISOString(timeInstance.offsetToDate(stopTimeOffset)); - CT::string resTime = "PT1S"; + CT::string resTime = "PT10M"; LayerMetadataDim dim = { .serviceName = "time", .cdfName = "time", .units = "ISO8601", .values = startTime + "/" + stopTime + "/" + resTime, .defaultValue = stopTime, .hasMultipleValues = true, .hidden = false}; metadataLayer->layerMetadata.dimList.push_back(dim); diff --git a/adagucserverEC/LayerTypeLiveUpdate/LayerTypeLiveUpdate.h b/adagucserverEC/LayerTypeLiveUpdate/LayerTypeLiveUpdate.h index 7326981b1..d68376e8c 100644 --- a/adagucserverEC/LayerTypeLiveUpdate/LayerTypeLiveUpdate.h +++ b/adagucserverEC/LayerTypeLiveUpdate/LayerTypeLiveUpdate.h @@ -11,11 +11,21 @@ */ int layerTypeLiveUpdateConfigureDimensionsInDataSource(CDataSource *dataSource); +/** + * Renders a Live Update layer depending on the presence or absence of post-processors. + */ +int layerTypeLiveUpdateRender(CDataSource *dataSource, CServerParams *srvParam); + /** * Renders time into an image the based on the time dimension in the GetMap request. */ int layerTypeLiveUpdateRenderIntoDrawImage(CDrawImage *image, CServerParams *srvParam); +/** + * Renders an image algorithmically based on the time dimension in the GetMap request. + */ +int layerTypeLiveUpdateRenderIntoImageDataWriter(CDataSource *dataSource, CServerParams *srvParam); + /** * Configures an actual time range in a WMSLayer object. This is used for generating the Layer element in the WMS GetCapabilities file */ diff --git a/adagucserverEC/utils/LayerMetadataStore.cpp b/adagucserverEC/utils/LayerMetadataStore.cpp index f707e2f07..85df254e8 100644 --- a/adagucserverEC/utils/LayerMetadataStore.cpp +++ b/adagucserverEC/utils/LayerMetadataStore.cpp @@ -289,7 +289,8 @@ int storeLayerStyleListIntoMetadataDb(MetadataLayer *metadataLayer) { } int loadLayerStyleListFromMetadataDb(MetadataLayer *metadataLayer) { - if (metadataLayer->dataSource->dLayerType == CConfigReaderLayerTypeCascaded || metadataLayer->dataSource->dLayerType == CConfigReaderLayerTypeLiveUpdate) { + if (metadataLayer->dataSource->dLayerType == CConfigReaderLayerTypeCascaded || + (metadataLayer->dataSource->dLayerType == CConfigReaderLayerTypeLiveUpdate) && metadataLayer->dataSource->cfgLayer->DataPostProc.empty()) { return 0; } if (!metadataLayer->readFromDb) { diff --git a/adagucserverEC/utils/XMLGenUtils.cpp b/adagucserverEC/utils/XMLGenUtils.cpp index 5b29ce9cb..5a7a31b1a 100644 --- a/adagucserverEC/utils/XMLGenUtils.cpp +++ b/adagucserverEC/utils/XMLGenUtils.cpp @@ -650,7 +650,9 @@ int getProjectionInformationForLayer(MetadataLayer *metadataLayer) { } int getStylesForLayer(MetadataLayer *metadataLayer) { - if (metadataLayer->dataSource->dLayerType == CConfigReaderLayerTypeCascaded || metadataLayer->dataSource->dLayerType == CConfigReaderLayerTypeLiveUpdate) { + if (metadataLayer->dataSource->dLayerType == CConfigReaderLayerTypeCascaded || + (metadataLayer->dataSource->dLayerType == CConfigReaderLayerTypeLiveUpdate && metadataLayer->dataSource->cfgLayer->DataPostProc.empty())) { + // Ignore styling in default case of the demo liveupdate layer, but not if there are data postprocessors return 0; } @@ -750,7 +752,8 @@ int getFileNameForLayer(MetadataLayer *metadataLayer) { } CServerParams *srvParam = metadataLayer->dataSource->srvParams; - if (metadataLayer->dataSource->dLayerType == CConfigReaderLayerTypeDataBase || metadataLayer->dataSource->dLayerType == CConfigReaderLayerTypeStyled) { + if (metadataLayer->dataSource->dLayerType == CConfigReaderLayerTypeDataBase || metadataLayer->dataSource->dLayerType == CConfigReaderLayerTypeStyled || + metadataLayer->dataSource->dLayerType == CConfigReaderLayerTypeLiveUpdate) { if (metadataLayer->dataSource->cfgLayer->Dimension.size() == 0) { metadataLayer->fileName.copy(metadataLayer->dataSource->cfgLayer->FilePath[0]->value.c_str()); if (CAutoConfigure::autoConfigureDimensions(metadataLayer->dataSource) != 0) { diff --git a/data/config/datasets/adaguc.tests.solarterminator.xml b/data/config/datasets/adaguc.tests.solarterminator.xml new file mode 100644 index 000000000..edf9028b5 --- /dev/null +++ b/data/config/datasets/adaguc.tests.solarterminator.xml @@ -0,0 +1,27 @@ + + + + + + {ADAGUC_PATH}data/datasets/solt.nc + + solarterminator + soltstyle + + + + \ No newline at end of file diff --git a/data/datasets/solt.nc b/data/datasets/solt.nc new file mode 100644 index 000000000..30563b6b3 Binary files /dev/null and b/data/datasets/solt.nc differ diff --git a/doc/configuration/2024-10-11-solt-layer.png b/doc/configuration/2024-10-11-solt-layer.png new file mode 100644 index 000000000..d771eeb93 Binary files /dev/null and b/doc/configuration/2024-10-11-solt-layer.png differ diff --git a/doc/configuration/Layer.md b/doc/configuration/Layer.md index 9b329ca98..7a7855b23 100644 --- a/doc/configuration/Layer.md +++ b/doc/configuration/Layer.md @@ -71,3 +71,38 @@ This layer type displays a GetMap image with the current time per second for the ``` + +If you configure the ```solarterminator``` data postprocessor, the liveupdate layer will display a GetMap image showing the areas where it is day, night, and different twilight levels (values from 0 to 4). The solar terminator algorithm is based on the calculations presented in the book *Astronomical Algorithms* (1991) by Jean Meeus. + +More information on how to configure this type of layer can be found [here](../tutorials/Configure_solar_terminator.md). + + +```xml + + + + + + {ADAGUC_PATH}data/datasets/solt.nc + + solarterminator + soltstyle + + +``` + + \ No newline at end of file diff --git a/doc/tutorials/Configure_solar_terminator.md b/doc/tutorials/Configure_solar_terminator.md new file mode 100644 index 000000000..94dfd34e5 --- /dev/null +++ b/doc/tutorials/Configure_solar_terminator.md @@ -0,0 +1,66 @@ +# Configure a Solar Terminator layer using AdagucServer + +- [Back to readme](./Readme.md) + +## Prerequisites + +Make sure adaguc-server is running, follow the instructions at [Starting the adaguc-server with docker](../Running.md) + +## Step 1: Copy a dummy file with timesteps into the adaguc-data folder + +Copy the file solt.nc into the adaguc-data folder: + + +``` +cp ${ADAGUC_PATH}/data/datasets/solt.nc ${ADAGUC_DATA_DIR} +``` +This file is available in the adaguc-server repository with location `data/datasets/solt.nc` + +## Step 2: Configure a dataset for this datafile + +Create the following file at the filepath `$ADAGUC_DATASET_DIR/solt.xml` + +```xml + + + + + + {ADAGUC_PATH}data/datasets/solt.nc + + solarterminator + soltstyle + + + + + +``` + +Feel free to customize the style to your preference, such as adding a fill color with transparency. However, keep in mind that the solar terminator calculation returns fixed categories rather than precise solar zenith angle values. + +## Step 3: Scan the new data + +``` +docker exec -i -t my-adaguc-server /adaguc/adaguc-server-updatedatasets.sh solt +``` + +## Step 4: Check if the layer works +You can use the following URL to try out your local Solar Terminator layer. + +Visit: +- https://adaguc.knmi.nl/adaguc-viewer/index.html?autowms=https://yourhostname/autowms#addlayer('https://yourhostname/adagucserver?dataset=solt&','solarterminator') diff --git a/doc/tutorials/Readme.md b/doc/tutorials/Readme.md index 47ddc60f3..b3cb4eec1 100644 --- a/doc/tutorials/Readme.md +++ b/doc/tutorials/Readme.md @@ -5,6 +5,7 @@ Tutorials - [Create a WMS service on a series of KNMI HDF5 files obtained from the KNMI data platform](Create_WMS_on_series_KNMI_HDF5_radar_files.md) - [Configure EDR timeseries service](Configure_EDR_service.md) +- [Configure a solar terminator layer](Configure_solar_terminator.md)