diff --git a/core/avhrrsatellite.cpp b/core/avhrrsatellite.cpp index 77529ce..c0c8a43 100644 --- a/core/avhrrsatellite.cpp +++ b/core/avhrrsatellite.cpp @@ -332,7 +332,7 @@ void AVHRRSatellite::AddSegmentsToList(QFileInfoList fileinfolist) } else delete segviirsdnbnoaa20; - } else if (fileInfo.fileName().mid( 0, 12) == "S3A_OL_1_EFR") // && fileInfo.fileName().mid( 100, 3) == "tar") // S3A EFR + } else if (fileInfo.fileName().mid( 0, 12) == "S3A_OL_1_EFR" || fileInfo.fileName().mid( 0, 12) == "S3B_OL_1_EFR") // S3A/B EFR { //S3A_OL_1_EFR____20161026T121318_20161026T121318_20161026T163853_0000_010_166______MAR_O_NR_002.SEN3.tar //0123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012 @@ -347,7 +347,7 @@ void AVHRRSatellite::AddSegmentsToList(QFileInfoList fileinfolist) } else delete segolciefr; - } else if (fileInfo.fileName().mid( 0, 12) == "S3A_OL_1_ERR") // && fileInfo.fileName().mid( 100, 3) == "tar") // S3A ERR + } else if (fileInfo.fileName().mid( 0, 12) == "S3A_OL_1_ERR" || fileInfo.fileName().mid( 0, 12) == "S3B_OL_1_ERR") // S3A/B ERR { //S3A_OL_1_ERR____20161026T121318_20161026T121318_20161026T163853_0000_010_166______MAR_O_NR_002.SEN3.tar //0123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012 @@ -362,7 +362,7 @@ void AVHRRSatellite::AddSegmentsToList(QFileInfoList fileinfolist) } else delete segolcierr; - } else if (fileInfo.fileName().mid( 0, 12) == "S3A_SL_1_RBT") + } else if (fileInfo.fileName().mid( 0, 12) == "S3A_SL_1_RBT" || fileInfo.fileName().mid( 0, 12) == "S3B_SL_1_RBT") { //S3A_SL_1_RBT____20170212T114405_20170212T114705_20170212T135851_0179_014_180_1800_SVL_O_NR_002.zip //0123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012 @@ -779,8 +779,8 @@ void AVHRRSatellite::ReadDirectories(QDate seldate, int hoursbefore) fileinfolist = map.values(); -// for(int i = 0; i < fileinfolist.count(); i++) -// qDebug() << "map values = " << fileinfolist.at(i).absoluteFilePath(); + for(int i = 0; i < fileinfolist.count(); i++) + qDebug() << "map values = " << fileinfolist.at(i).absoluteFilePath(); emit signalResetProgressbar(fileinfolist.size(), (*its)); @@ -1140,7 +1140,7 @@ void AVHRRSatellite::InsertToMap(QFileInfoList fileinfolist, QMap 0); m_pBuffer->clear(); @@ -77,6 +81,9 @@ void DatahubAccessManager::DownloadXML(QDate selectdate, eDatahub hub) QUrl url; qDebug() << "start DownloadXML"; + qDebug()<<"SSL version use for build: "< cnt_spectrum(nbr_spectrum); - widget->clear(); QMap > >::const_iterator citdate = map.constBegin(); @@ -346,7 +345,6 @@ void FormGeostationary::slotCreateGeoImage(QString type, QVector spectr QString tex; - SegmentListGeostationary *sl; for(int i = 0; i < opts.geosatellites.count(); i++) @@ -361,16 +359,9 @@ void FormGeostationary::slotCreateGeoImage(QString type, QVector spectr } } -// int geoindex = formtoolbox->getGeoIndex(); int geoindex = sl->getGeoSatelliteIndex(); - formtoolbox->createFilenamestring(opts.geosatellites.at(geoindex).shortname, tex, spectrumvector); - - sl = setActiveSegmentList(geoindex); sl->ResetSegments(); - - qDebug() << "FormGeostationary::slotCreateGeoImage(eGeoSatellite, QString, QVector, QVector, int) geoindex " << geoindex; - imageptrs->ResetPtrImage(); qDebug() << QString(" CreateGeoImage ; kind of image = %1 spectrumvector = %2 ; %3 ; %4").arg(sl->getKindofImage()).arg(spectrumvector.at(0)).arg(spectrumvector.at(1)).arg(spectrumvector.at(2)); @@ -383,20 +374,20 @@ void FormGeostationary::slotCreateGeoImage(QString type, QVector spectr if (type == "HRV" || type == "HRV Color") { - if (sl->areatype == 1) - { - if(opts.geosatellites.at(geoindex).imageheighthrv1 > 0 && opts.geosatellites.at(geoindex).imagewidthhrv1) - imageptrs->InitializeImageGeostationary(opts.geosatellites.at(geoindex).imagewidthhrv1, opts.geosatellites.at(geoindex).imageheighthrv1); - else - return; - } + if (sl->areatype == 1) + { + if(opts.geosatellites.at(geoindex).imageheighthrv1 > 0 && opts.geosatellites.at(geoindex).imagewidthhrv1) + imageptrs->InitializeImageGeostationary(opts.geosatellites.at(geoindex).imagewidthhrv1, opts.geosatellites.at(geoindex).imageheighthrv1); else - { - if(opts.geosatellites.at(geoindex).imageheighthrv0 > 0 && opts.geosatellites.at(geoindex).imagewidthhrv0 > 0) - imageptrs->InitializeImageGeostationary(opts.geosatellites.at(geoindex).imagewidthhrv0, opts.geosatellites.at(geoindex).imageheighthrv0); - else - return; - } + return; + } + else + { + if(opts.geosatellites.at(geoindex).imageheighthrv0 > 0 && opts.geosatellites.at(geoindex).imagewidthhrv0 > 0) + imageptrs->InitializeImageGeostationary(opts.geosatellites.at(geoindex).imagewidthhrv0, opts.geosatellites.at(geoindex).imageheighthrv0); + else + return; + } } else { @@ -406,16 +397,14 @@ void FormGeostationary::slotCreateGeoImage(QString type, QVector spectr return; } - formimage->displayImage(IMAGE_GEOSTATIONARY); - - //formimage->adjustPicSize(true); + formimage->adjustPicSize(false); qDebug() << QString("FormGeostationary::CreateGeoImage kind = %1 areatype = %2").arg(type).arg(sl->areatype); - if(opts.geosatellites.at(sl->getGeoSatelliteIndex()).protocol == "HDF" ) + if(opts.geosatellites.at(geoindex).protocol == "HDF" ) CreateGeoImageHDF(sl, type, tex, spectrumvector, inversevector); - else if(opts.geosatellites.at(sl->getGeoSatelliteIndex()).protocol == "netCDF") + else if(opts.geosatellites.at(geoindex).protocol == "netCDF") CreateGeoImagenetCDF(sl, type, tex, spectrumvector, inversevector, histogrammethod, pseudocolor); else CreateGeoImageXRIT(sl, type, tex, spectrumvector, inversevector, histogrammethod); @@ -547,7 +536,6 @@ void FormGeostationary::CreateGeoImageXRIT(SegmentListGeostationary *sl, QString { if(opts.geosatellites.at(geoindex).prologfile) { - fa = (type == "VIS_IR" || type == "VIS_IR Color" ? faVIS_IR : faHRV); // Read prologue @@ -567,8 +555,6 @@ void FormGeostationary::CreateGeoImageXRIT(SegmentListGeostationary *sl, QString qDebug() << QString("Error : runtime error in reading prologue file : %1").arg(run.what()); } } - - } if(opts.geosatellites.at(geoindex).epilogfile) @@ -620,7 +606,6 @@ void FormGeostationary::CreateGeoImageXRIT(SegmentListGeostationary *sl, QString sl->UpperSouthLineActual = 0; sl->UpperWestColumnActual = 0; sl->UpperNorthLineActual = 0; - } } } @@ -914,633 +899,28 @@ int FormGeostationary::getTabWidgetGeoIndex() void FormGeostationary::slotCreateRGBrecipe(int recipe) { - qDebug() << "CreateRGBrecipe = " << recipe; - CreateRGBrecipeImage(recipe); -} - -void FormGeostationary::CreateRGBrecipeImage(int recipe) -{ + QRgb *row_col; - QString tex; - QStringList redbandlist = imageptrs->rgbrecipes[recipe].Colorvector.at(0).channels; - QStringList greenbandlist = imageptrs->rgbrecipes[recipe].Colorvector.at(1).channels; - QStringList bluebandlist = imageptrs->rgbrecipes[recipe].Colorvector.at(2).channels; - QStringList uniquebandlist; - QList uniqueunitlist; - seviri_units redunits = imageptrs->rgbrecipes[recipe].Colorvector.at(0).units; - seviri_units greenunits = imageptrs->rgbrecipes[recipe].Colorvector.at(1).units; - seviri_units blueunits = imageptrs->rgbrecipes[recipe].Colorvector.at(2).units; - QList uniquechannelnbrlist; - - for(int i = 0; i < redbandlist.length(); i++) - { - if(!uniquebandlist.contains(redbandlist.at(i))) - { - uniquebandlist.append(redbandlist.at(i)); - uniqueunitlist.append(redunits); - uniquechannelnbrlist.append(imageptrs->rgbrecipes[recipe].Colorvector.at(0).spectral_channel_nbr.at(i)); - } - } - for(int i = 0; i < greenbandlist.length(); i++) - { - if(!uniquebandlist.contains(greenbandlist.at(i))) - { - uniquebandlist.append(greenbandlist.at(i)); - uniqueunitlist.append(greenunits); - uniquechannelnbrlist.append(imageptrs->rgbrecipes[recipe].Colorvector.at(1).spectral_channel_nbr.at(i)); - } - } - for(int i = 0; i < bluebandlist.length(); i++) - { - if(!uniquebandlist.contains(bluebandlist.at(i))) - { - uniquebandlist.append(bluebandlist.at(i)); - uniqueunitlist.append(blueunits); - uniquechannelnbrlist.append(imageptrs->rgbrecipes[recipe].Colorvector.at(2).spectral_channel_nbr.at(i)); - } - } + qDebug() << "CreateRGBrecipe = " << recipe; - QString recipename = imageptrs->rgbrecipes[recipe].Name; - qDebug() << "recipename = " << recipename; - qDebug() << "uniquebandlist contains = " << uniquebandlist.count(); - for(int i = 0; i < uniquebandlist.length(); i++) - { - qDebug() << "channels at " << uniquechannelnbrlist.at(i) << " = " << uniquebandlist.at(i); - } + //CreateRGBrecipeImageThread(recipe); SegmentListGeostationary *sl; + QString tex; for(int i = 0; i < opts.geosatellites.count(); i++) { QList treewidgetselected = geotreewidgetlist.at(i)->selectedItems(); if(treewidgetselected.count() > 0) { - sl = setActiveSegmentList(i); QTreeWidgetItem *it = treewidgetselected.at(0); tex = it->text(0); + sl = setActiveSegmentList(i); break; } } - qDebug() << "tex = " << tex << " imagepath = " << sl->getImagePath(); - //tex = "2018-02-13 14:45" - - int geoindex = sl->getGeoSatelliteIndex(); - sl = setActiveSegmentList(geoindex); - sl->ResetSegments(); - imageptrs->ResetPtrImage(); - - QString directory = sl->getImagePath(); - QString productid1 = opts.geosatellites.at(geoindex).searchstring.mid(6, 4); - QString productid2 = "IR_016"; - QString timing = tex.mid(0, 4) + tex.mid(5, 2) + tex.mid(8, 2) + tex.mid(13, 2) + tex.mid(16, 2); - - qDebug() << "directory = " << directory; - qDebug() << "productid1 = " << productid1; - qDebug() << "productid2 = " << productid2; - qDebug() << "timing = " << timing; - - int totalpixels = 3712 * 3712; - - double *time = new double[totalpixels]; - float *lat = new float[totalpixels]; - float *lon = new float[totalpixels]; - float *sza = new float[totalpixels]; -// float *saa = new float[totalpixels]; -// float *vza = new float[totalpixels]; -// float *vaa = new float[totalpixels]; - - QList result; - result.append(new float[totalpixels]); - result.append(new float[totalpixels]); - result.append(new float[totalpixels]); - float resultmax[3]; - float resultmin[3]; - - for(int i = 0; i < 3; i++) - { - resultmin[i] = std::numeric_limits::max(); - resultmax[i] = std::numeric_limits::min(); - } - - snu_init_array_d(time, totalpixels, FILL_VALUE_F); - snu_init_array_f(lat, totalpixels, FILL_VALUE_F); - snu_init_array_f(lon, totalpixels, FILL_VALUE_F); - snu_init_array_f(sza, totalpixels, FILL_VALUE_F); -// snu_init_array_f(saa, totalpixels, FILL_VALUE_F); -// snu_init_array_f(vza, totalpixels, FILL_VALUE_F); -// snu_init_array_f(vaa, totalpixels, FILL_VALUE_F); - snu_init_array_f(result[0], totalpixels, FILL_VALUE_F); - snu_init_array_f(result[1], totalpixels, FILL_VALUE_F); - snu_init_array_f(result[2], totalpixels, FILL_VALUE_F); - -// snu_init_array_f(d2->data2, d->image.n_bands * length, d2->fill_value); - - typedef struct { - int spectral_channel_nbr; - QString spectrum; - float min; - float max; - float *data; - seviri_units units; - } bandstorage; - - QList bands; - - for(int i = 0; i < uniquebandlist.length(); i++) - { - bandstorage newband; - newband.spectral_channel_nbr = uniquechannelnbrlist.at(i); - newband.data = new float[totalpixels]; - snu_init_array_f(newband.data, totalpixels, FILL_VALUE_F); - newband.min = 0; newband.max = 0; - newband.spectrum = uniquebandlist.at(i); - newband.units = uniqueunitlist.at(i); - bands.append(newband); - } - - qDebug() << "number of bands = " << bands.length(); - - - MsgFileAccess fa(directory, "H", productid1, productid2, timing); - MsgDataAccess da; - MSG_data prodata; - MSG_data epidata; - MSG_header header; - - - da.scan(fa, prodata, epidata, header); - - long navloff = header.image_navigation->LOFF; - long navcoff = header.image_navigation->COFF; - long navlfac = header.image_navigation->LFAC; - long navcfac = header.image_navigation->CFAC; - - qDebug() << QString("image navigation LOFF = %1").arg(navloff); //, 16, 'f', 2); - qDebug() << QString("image navigation COFF = %1").arg(navcoff); //, 16, 'f', 2); - qDebug() << QString("image navigation LFAC = %1").arg(navlfac); //, 16, 'f', 2); - qDebug() << QString("image navigation CFAC = %1").arg(navcfac); //, 16, 'f', 2); - - qDebug() << QString("column scaling factor = %1").arg(header.image_navigation->column_scaling_factor); //, 16, 'f', 2); - qDebug() << QString("line scaling factor = %1").arg(header.image_navigation->line_scaling_factor); //, 16, 'f', 2); - qDebug() << QString("column offset = %1").arg(header.image_navigation->column_offset); //, 16, 'f', 2); - qDebug() << QString("line offset = %1").arg(header.image_navigation->line_offset); //, 16, 'f', 2); - - - //cout << epidata; - - qDebug()<< epidata.epilogue->product_stats.ActualScanningSummary.ForwardScanStart.get_timestring().c_str(); - double jtime_start = epidata.epilogue->product_stats.ActualScanningSummary.ForwardScanStart.get_jtime(); - double jtime_end = epidata.epilogue->product_stats.ActualScanningSummary.ForwardScanEnd.get_jtime(); - - double jtime2; - double jtime = (jtime_start + jtime_end) / 2.; - - struct tm cdate; - Calendar_Date(jtime, &cdate); - double day_of_year = jtime - QSgp4Date::DateToJD(cdate.tm_year, 1, 0, true); - - qDebug() << "day_of_year = " << day_of_year; - - double jtime_start2; - double jtime_end2; - - /*------------------------------------------------------------------------- - * Compute the satellite position vector in Cartesian coordinates (km). - *-----------------------------------------------------------------------*/ - int i; - for (i = 0; i < 100; ++i) { - jtime_start2 = prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].StartTime.get_jtime(); - jtime_end2 = prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].EndTime.get_jtime(); - if (jtime >= jtime_start2 && jtime <= jtime_end2) - break; - } - - if (i == 100) { - fprintf(stderr, "ERROR: Image time is out of range of supplied orbit " - "polynomials\n"); - return; - } - - double t, X, Y, Z; - t = (jtime - (jtime_start2 + jtime_end2) / 2.) / ((jtime_end2 - jtime_start2) / 2.); - - X = prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].X[0] + - prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].X[1] * t; - Y = prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].Y[0] + - prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].Y[1] * t; - Z = prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].Z[0] + - prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].Z[1] * t; - - qDebug() << "i = " << i; - qDebug() << "X = " << X; - qDebug() << "Y = " << Y; - qDebug() << "Z = " << Z; - - /*------------------------------------------------------------------------- - * Compute latitude and longitude and solar and sensor zenith and azimuth - * angles. - *-----------------------------------------------------------------------*/ - - double lon0 = prodata.prologue->image_description.ProjectionDescription.LongitudeOfSSP; - qDebug() << "longitude SSP = " << lon0; - - double mu0; - double theta0; - double phi0; - - long countF = 0; - - quint16 seglines = header.image_structure->number_of_lines; - quint16 columns = header.image_structure->number_of_columns; - - qDebug() << "Start calculations"; - unsigned int i_image; - - - for(int i = 0; i < uniquebandlist.length(); i++) - { - MsgFileAccess fac(fa, uniquebandlist.at(i)); - da.scan(fac, header); - bands[i].spectral_channel_nbr = header.segment_id->spectral_channel_id; - MSG_SAMPLE compmin = 0xffff, compmax = 0; - for (size_t j = 0; j < da.segnames.size(); ++j) - { - cout << "Segment " << j << ": "; - MSG_data* d = da.segment(j); - MSG_SAMPLE min = 0xffff, max = 0; - if (!d) - { - cout << "missing.\n"; - continue; - } else { - for (size_t k = 0; k < da.npixperseg; ++k) - { - if (d->image->data[k] < min) min = d->image->data[k]; - if (d->image->data[k] > max) max = d->image->data[k]; - bands[i].data[j*da.npixperseg + k] = d->image->data[k]; - } - cout << "min " << min << " max " << max << endl; - } - if (min < compmin) compmin = min; - if (max > compmax) compmax = max; - } - cout << "compmin = " << compmin << " compmax = " << compmax << endl; - bands[i].min = (float)compmin; - bands[i].max = (float)compmax; - } - - - for (int i = 0; i < 3712; ++i) - { - jtime2 = jtime_start + (double) i / (double) (3712 - 1) * (jtime_end - jtime_start); - for (int j = 0; j < 3712; ++j) - { - i_image = i * 3712 + j; - - snu_line_column_to_lat_lon(i, j, &lat[i_image], &lon[i_image], lon0, &nav_scaling_factors_vir); - - if (lat[i_image] != FILL_VALUE_F && lon[i_image] != FILL_VALUE_F) - { - time[i_image] = jtime2; - - snu_solar_params2(jtime2, lat[i_image] * D2R, lon[i_image] * D2R, &mu0, &theta0, &phi0, NULL); - sza[i_image] = theta0 * R2D; -// saa[i_image] = phi0 * R2D; - -// saa[i_image] = saa[i_image] + 180.; -// if (saa[i_image] > 360.) -// saa[i_image] = saa[i_image] - 360.; - -// snu_vza_and_vaa(lat[i_image], lon[i_image], 0., X, Y, Z, &vza[i_image], &vaa[i_image]); - -// vaa[i_image] = vaa[i_image] + 180.; -// if (vaa[i_image] > 360.) -// vaa[i_image] = vaa[i_image] - 360.; - } - else - { - countF++; - for(int k = 0; k < bands.length(); k++) - bands[k].data[i_image] = FILL_VALUE_F; - } - } - } - qDebug() << "End calculations countF = " << countF; - - - for(int i = 0; i < bands.length(); i++) - { - qDebug() << "spectral channel = " << bands[i].spectral_channel_nbr; - qDebug() << "spectrum = " << bands[i].spectrum; - qDebug() << "min = " << bands[i].min << " max = " << bands[i].max; - qDebug() << "units = " << (int)bands[i].units; - } - - /*------------------------------------------------------------------------- - * Compute radiance for the bands requested. - * - * Ref: PDF_TEN_05105_MSG_IMG_DATA, Page 26 - *-----------------------------------------------------------------------*/ - - double slope, offset; - bool toint; - - - for(int i = 0; i < bands.length(); i++) - { - if(bands[i].units == SEVIRI_UNIT_RAD) // mW*pow(m,-2)*pow(sr,-1)*pow(pow(cm,-1)), -1) - { - prodata.prologue->radiometric_proc.get_slope_offset(bands[i].spectral_channel_nbr, slope, offset, toint); - bands[i].min = std::numeric_limits::max(); - bands[i].max = std::numeric_limits::min(); - for (int j = 0; j < 3712; ++j) - { - for (int k = 0; k < 3712; ++k) - { - i_image = j * 3712 + k; - - if (bands[i].data[i_image] != FILL_VALUE_US && bands[i].data[i_image] > 0) - { - bands[i].data[i_image] = bands[i].data[i_image] * slope + offset; - if(bands[i].data[i_image] < bands[i].min) bands[i].min = bands[i].data[i_image]; - if(bands[i].data[i_image] > bands[i].max) bands[i].max = bands[i].data[i_image]; - } - } - } - } - } - - /*------------------------------------------------------------------------- - * Compute reflectance or bidirectional reflectance factor (BRF) for the - * bands requested. - * - * Ref: PDF_MSG_SEVIRI_RAD2REFL, Page 8 - *-----------------------------------------------------------------------*/ - double dd = 1. / sqrt(snu_solar_distance_factor2(day_of_year)); - - double a = PI * dd * dd; - - int satid = (int)header.segment_id->spacecraft_id - 321 ; - - for (int i = 0; i < bands.length(); ++i) - { - if (bands[i].units == SEVIRI_UNIT_REF || bands[i].units == SEVIRI_UNIT_BRF) - { - prodata.prologue->radiometric_proc.get_slope_offset(bands[i].spectral_channel_nbr, slope, offset, toint); - bands[i].min = std::numeric_limits::max(); - bands[i].max = std::numeric_limits::min(); - - double b = a / band_solar_irradiance[satid][bands[i].spectral_channel_nbr - 1]; - - for (int j = 0; j < 3712; ++j) - { - for (int k = 0; k < 3712; ++k) - { - i_image = j * 3712 + k; - - if (bands[i].units == SEVIRI_UNIT_BRF && bands[i].data[i_image] != FILL_VALUE_US && bands[i].data[i_image] > 0 && - sza[i_image] >= 0. && sza[i_image] < 80.) - { - double R = bands[i].data[i_image] * slope + offset; - - bands[i].data[i_image] = b * R; - - if (bands[i].units == SEVIRI_UNIT_BRF) - bands[i].data[i_image] /= cos(sza[i_image] * D2R); - - if(bands[i].data[i_image] < bands[i].min) bands[i].min = bands[i].data[i_image]; - if(bands[i].data[i_image] > bands[i].max) bands[i].max = bands[i].data[i_image]; - } - else if (bands[i].units == SEVIRI_UNIT_REF && bands[i].data[i_image] != FILL_VALUE_US && bands[i].data[i_image] > 0) - { - double R = bands[i].data[i_image] * slope + offset; - - bands[i].data[i_image] = b * R; - - if(bands[i].data[i_image] < bands[i].min) bands[i].min = bands[i].data[i_image]; - if(bands[i].data[i_image] > bands[i].max) bands[i].max = bands[i].data[i_image]; - } - } - } - } - } - - const double c1 = 1.19104e-5; - const double c2 = 1.43877; - - /*------------------------------------------------------------------------- - * Compute brightness temperature for the bands requested. - * - * Ref: PDF_TEN_05105_MSG_IMG_DATA, Page 26 - * Ref: The Conversion from Effective Radiances - * to Equivalent Brightness Temperatures (EUM/MET/TEN/11/0569) - *-----------------------------------------------------------------------*/ - for (int i = 0; i < bands.length(); ++i) - { - if (bands[i].units == SEVIRI_UNIT_BT) - { - prodata.prologue->radiometric_proc.get_slope_offset(bands[i].spectral_channel_nbr, slope, offset, toint); - bands[i].min = std::numeric_limits::max(); - bands[i].max = std::numeric_limits::min(); - - qDebug() << bands[i].spectral_channel_nbr << " " << bands[i].spectrum << " slope = " << slope << " offset = " << offset; -/* - nu = 1.e4 / channel_center_wavelength[d->image.band_ids[i] - 1]; -*/ - double nu = bt_nu_c[satid][bands[i].spectral_channel_nbr - 1]; - double nu3 = nu * nu * nu; - - double a = bt_A[satid][bands[i].spectral_channel_nbr - 1]; - double b = bt_B[satid][bands[i].spectral_channel_nbr - 1]; - - for (int j = 0; j < 3712; ++j) - { - for (int k = 0; k < 3712; ++k) - { - i_image = j * 3712 + k; - - if (bands[i].data[i_image] != FILL_VALUE_F && bands[i].data[i_image] > 0) - { - double L = bands[i].data[i_image] * slope + offset; - bands[i].data[i_image] = (c2 * nu / log(1. + nu3 * c1 / L) - b) / a; - if(bands[i].data[i_image] < bands[i].min) bands[i].min = bands[i].data[i_image]; - if(bands[i].data[i_image] > bands[i].max) bands[i].max = bands[i].data[i_image]; - } - } - } - } - } - - QRgb *row_col; - float red, green, blue; - -// for (int line = 0; line < 3712; line++) -// { -// row_col = (QRgb*)imageptrs->ptrimageGeostationary->scanLine(3711 - line); -// for (int pixelx = 0; pixelx < 3712; pixelx++) -// { -// i_image = line * 3712 + pixelx; -// if(bands[0].data[i_image] != FILL_VALUE_F ) -// { -// red = 255.0 * pow((bands[0].data[i_image] - bands[0].min) / (bands[0].max - bands[0].min), 1.0/imageptrs->rgbrecipes[recipe].Colorvector.at(0).gamma); -// green = 255.0 * pow((bands[1].data[i_image] - bands[1].min) / (bands[1].max - bands[1].min), 1.0/imageptrs->rgbrecipes[recipe].Colorvector.at(1).gamma); -// blue = 255.0 * pow((bands[2].data[i_image] - bands[2].min) / (bands[2].max - bands[2].min), 1.0/imageptrs->rgbrecipes[recipe].Colorvector.at(2).gamma); -// } -// else -// { -// red = 0.0; -// green = 0.0; -// blue = 0.0; -// } -// row_col[3711 - pixelx] = qRgb((int)red, (int)green, (int)blue); -// } -// } - - countF = 0; - for(int colorindex = 0; colorindex < 3; colorindex++) - { - for(int i = 0; i < imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).channels.length(); i++) - { - for(int j = 0; j < bands.length(); j++) - { - if(bands[j].spectral_channel_nbr == imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).spectral_channel_nbr.at(i)) - { - qDebug() << colorindex << " " << imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).channels.at(i) << " " << - imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).spectral_channel_nbr.at(i) << " " << - imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).subtract.at(i); - for (int line = 0; line < 3712; line++) - { - row_col = (QRgb*)imageptrs->ptrimageGeostationary->scanLine(3711 - line); - for (int pixelx = 0; pixelx < 3712; pixelx++) - { - i_image = line * 3712 + pixelx; - if(bands[j].data[i_image] != FILL_VALUE_F ) - { - if(imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).subtract.at(i)) - { - if(result[colorindex][i_image] == FILL_VALUE_F) - result[colorindex][i_image] = - bands[j].data[i_image]; - else - result[colorindex][i_image] -= bands[j].data[i_image]; - } - else - { - if(result[colorindex][i_image] == FILL_VALUE_F) - result[colorindex][i_image] = bands[j].data[i_image]; - else - result[colorindex][i_image] += bands[j].data[i_image]; - } - if(imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).dimension == "K") - { - if(result[colorindex][i_image] > imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).rangeto ) - result[colorindex][i_image] = imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).rangeto; - if(result[colorindex][i_image] < imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).rangefrom ) - result[colorindex][i_image] = imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).rangefrom; - } - if(result[colorindex][i_image] < resultmin[colorindex]) resultmin[colorindex] = result[colorindex][i_image]; - if(result[colorindex][i_image] > resultmax[colorindex]) resultmax[colorindex] = result[colorindex][i_image]; - - } - else - { - countF++; - result[colorindex][i_image] = FILL_VALUE_F; - } - } - } - } - } - } - } - - qDebug() << "countF = " << countF; - - for(int colorindex = 0; colorindex < 3; colorindex++ ) - qDebug() << QString("%1 resultmin = %2 resultmax = %3").arg(colorindex).arg(resultmin[colorindex]).arg(resultmax[colorindex]); - - for (int line = 0; line < 3712; line++) - { - row_col = (QRgb*)imageptrs->ptrimageGeostationary->scanLine(3711 - line); - for (int pixelx = 0; pixelx < 3712; pixelx++) - { - i_image = line * 3712 + pixelx; - - if(result[0][i_image] != FILL_VALUE_F) // && result[1][i_image] != FILL_VALUE_F && result[2][i_image] != FILL_VALUE_F ) - { - red = 255.0 * pow((result[0][i_image] - resultmin[0]) / (resultmax[0] - resultmin[0]), 1.0/imageptrs->rgbrecipes[recipe].Colorvector.at(0).gamma); - green = 255.0 * pow((result[1][i_image] - resultmin[1]) / (resultmax[1] - resultmin[1]), 1.0/imageptrs->rgbrecipes[recipe].Colorvector.at(1).gamma); - blue = 255.0 * pow((result[2][i_image] - resultmin[2]) / (resultmax[2] - resultmin[2]), 1.0/imageptrs->rgbrecipes[recipe].Colorvector.at(2).gamma); - } - else - { - red = 0.0; - green = 0.0; - blue = 0.0; - } - row_col[3711 - pixelx] = qRgb((int)red, (int)green, (int)blue); - } - } - - - - - int i_pixel = 1856 * 3712 + 1856; - - qDebug() << QString("time = %1").arg(time[1856 * 3712 + 1856], 16, 'f', 5); - qDebug() << QString("lat = %1").arg(lat[1856 * 3712 + 1856], 16, 'f', 5); - qDebug() << QString("lon = %1").arg(lon[1856 * 3712 + 1856], 16, 'f', 5); - qDebug() << QString("sza = %1").arg(sza[1856 * 3712 + 1856], 16, 'f', 5); -// qDebug() << QString("saa = %1").arg(saa[1856 * 3712 + 1856], 16, 'f', 5); -// qDebug() << QString("vza = %1").arg(vza[1856 * 3712 + 1856], 16, 'f', 5); -// qDebug() << QString("vaa = %1").arg(vaa[1856 * 3712 + 1856], 16, 'f', 5); - -// printf("Julian Day Number: % .8e\n", time[i_pixel]); -// printf("latitude: % .8e\n", lat [i_pixel]); -// printf("longitude: % .8e\n", lon [i_pixel]); -// printf("solar zenith angle: % .8e\n", sza [i_pixel]); -// printf("solar azimuth angle: % .8e\n", saa [i_pixel]); -// printf("viewing zenith angle: % .8e\n", vza [i_pixel]); -// printf("viewing azimuth angle: % .8e\n", vaa [i_pixel]); - - qDebug() << "min float = " << std::numeric_limits::min(); - qDebug() << "max float = " << std::numeric_limits::max(); - -// qDebug() << QString("bands[0] = %1 [1856][1856] = %2 units = %3 min = %4 max = %5") -// .arg(bands[0].spectrum).arg(bands[0].data[1856 * 3712 + 1856], 9, 'f', 3).arg((int)bands[0].units).arg(bands[0].min).arg(bands[0].max); -// qDebug() << QString("bands[1] = %1 [1856][1856] = %2 units = %3 min = %4 max = %5") -// .arg(bands[1].spectrum).arg(bands[1].data[1856 * 3712 + 1856], 9, 'f', 3).arg((int)bands[1].units).arg(bands[1].min).arg(bands[1].max); -// qDebug() << QString("bands[2] = %1 [1856][1856] = %2 units = %3 min = %4 max = %5") -// .arg(bands[2].spectrum).arg(bands[2].data[1856 * 3712 + 1856], 9, 'f', 3).arg((int)bands[2].units).arg(bands[2].min).arg(bands[2].max); - - qDebug() << "CreateRGBrecipeImage(int recipe) Finished !!"; - -// cout << "Columns: " << da.columns << endl; -// cout << "Lines: " << da.lines << endl; -// cout << "Segments: " << da.segnames.size() << endl; -// cout << "Pixels per segment: " << da.npixperseg << endl; -// cout << "Scanlines per segment: " << da.seglines << endl; -// cout << "East is: " << (da.swapX ? "left" : "right") << endl; -// cout << "North is: " << (da.swapY ? "down" : "up") << endl; -// cout << "HRV: " << (da.hrv ? "yes" : "no") << endl; - - - for(int i = 0; i < uniquebandlist.length(); i++) - { - delete [] bands[i].data; - } - - delete [] time; - delete [] lat; - delete [] lon; - delete [] sza; -// delete [] saa; -// delete [] vza; -// delete [] vaa; - delete [] result[0]; - delete [] result[1]; - delete [] result[2]; - + sl->ComposeGeoRGBRecipe(recipe, tex); } diff --git a/core/formgeostationary.h b/core/formgeostationary.h index a1c055c..6112c16 100644 --- a/core/formgeostationary.h +++ b/core/formgeostationary.h @@ -10,6 +10,22 @@ #include "formtoolbox.h" #include "formimage.h" +//typedef struct { +// int listindex; +// int spectral_channel_nbr; +// QString directory; +// QString productid1; +// QString productid2; +// QString timing; +// double day_of_year; +// float min; +// float max; +// float *data; +// seviri_units units; +// double slope; +// double offset; +//} bandstorage; + class FormImage; namespace Ui { @@ -32,16 +48,17 @@ class FormGeostationary : public QWidget SegmentListGeostationary *setActiveSegmentList(int geoindex); SegmentListGeostationary *getActiveSegmentList(); int getTabWidgetGeoIndex(); + //void ComposeGeoRGBRecipe(bandstorage &bs); + ~FormGeostationary(); private: QStringList getGeostationarySegments(int geoindex, const QString imagetype, const QString filepath, QVector spectrumvector, QString filepattern); void PopulateTreeGeo(int geoindex); -// void PopulateTreeGeo1(int geoindex); void CreateGeoImageXRIT(SegmentListGeostationary *sl, QString type, QString tex, QVector spectrumvector, QVector inversevector, int histogrammethod); void CreateGeoImageHDF(SegmentListGeostationary *sl, QString type, QString tex, QVector spectrumvector, QVector inversevector); void CreateGeoImagenetCDF(SegmentListGeostationary *sl, QString type, QString tex, QVector spectrumvector, QVector inversevector, int histogrammethod, bool pseudocolor); - void CreateRGBrecipeImage(int recipe); + Ui::FormGeostationary *ui; AVHRRSatellite *segs; SatelliteList *sats; @@ -63,6 +80,7 @@ private slots: void geostationarysegmentschosen(int geoindex, QStringList ll); void setbuttonlabels(int geoindex, bool state); void enabletoolboxbuttons(bool); + }; #endif // FORMGEOSTATIONARY_H diff --git a/core/formimage.cpp b/core/formimage.cpp index ac4af62..beb4e4b 100644 --- a/core/formimage.cpp +++ b/core/formimage.cpp @@ -1376,13 +1376,17 @@ void FormImage::displayGeoImageInfo() void FormImage::displayGeoImageInformation(QString satname) { + QString spectrum; SegmentListGeostationary *slgeo = NULL; slgeo = segs->getActiveSegmentList(); QString type = slgeo->getKindofImage(); QVector spectrumvector = slgeo->getSpectrumVector(); - QString spectrum = ( type == "VIS_IR" ? spectrumvector.at(0) : ""); + if(spectrumvector.size() > 0) + spectrum = ( type == "VIS_IR" ? spectrumvector.at(0) : ""); + else + spectrum = ""; txtInfo = QString("" "Info" diff --git a/core/formmapcyl.ui b/core/formmapcyl.ui index f56e29a..2ab27d5 100644 --- a/core/formmapcyl.ui +++ b/core/formmapcyl.ui @@ -39,7 +39,7 @@ - 0 + 2 diff --git a/core/formtoolbox.cpp b/core/formtoolbox.cpp index d1296ca..4e9f543 100644 --- a/core/formtoolbox.cpp +++ b/core/formtoolbox.cpp @@ -132,7 +132,8 @@ FormToolbox::FormToolbox(QWidget *parent, FormImage *p_formimage, FormGeostation inversevector.append(false); formimage->channelshown = IMAGE_GEOSTATIONARY; - ui->tabWidget->setCurrentIndex(opts.currenttabwidget); + qDebug() << QString("Current tabwidget = %1").arg(opts.currenttabwidget); + //ui->tabWidget->setCurrentIndex(opts.currenttabwidget); ui->tabWidgetVIIRS->setCurrentIndex(0); ui->tabWidgetSentinel->setCurrentIndex(0); @@ -231,7 +232,6 @@ FormToolbox::FormToolbox(QWidget *parent, FormImage *p_formimage, FormGeostation qDebug() << QString("poi.strlComboGeo5.at(geoindex) = %1 ").arg(poi.strlComboGeo5.at(geoindex)); qDebug() << QString("poi.strlComboGeo6.at(geoindex) = %1 ").arg(poi.strlComboGeo6.at(geoindex)); - setPOIsettings(); setMConfigsettings(); setOLCIefrConfigsettings(); @@ -244,9 +244,6 @@ FormToolbox::FormToolbox(QWidget *parent, FormImage *p_formimage, FormGeostation qDebug() << QString("poi.strlComboGeo5.at(geoindex) = %1 ").arg(poi.strlComboGeo5.at(geoindex)); qDebug() << QString("poi.strlComboGeo6.at(geoindex) = %1 ").arg(poi.strlComboGeo6.at(geoindex)); - - - qDebug() << QString("Setting currenttoolbox = %1").arg(opts.currenttoolbox); ui->toolBox->setCurrentIndex(opts.currenttoolbox); // in projection tab LCC GVP or SG ui->comboPOI->setCurrentIndex(0); @@ -277,9 +274,6 @@ FormToolbox::FormToolbox(QWidget *parent, FormImage *p_formimage, FormGeostation ui->lblCentreBand->setText(QString("%1").arg(fval, 0, 'E', 2)); ui->lblTitleCentreBand->setText(QString("Centre Band from %1 to %2 [W/cm² sr]").arg(fval/pow(10, opts.dnbspbwindowsvalue), 0, 'E', 2).arg(fval*pow(10, opts.dnbspbwindowsvalue), 0, 'E', 2)); - - - ui->lblGeo1->setText("0.635"); ui->lblGeo2->setText("0.81"); ui->lblGeo3->setText("1.64"); @@ -517,8 +511,15 @@ bool FormToolbox::eventFilter(QObject *target, QEvent *event) } +void FormToolbox::setValuePrgBar(int val) +{ + ui->pbProgress->setValue(val); + //this->update(); +} + void FormToolbox::setValueProgressBar(int val) { + qDebug() << "FormToolbox::setValueProgressBar(int val) " << val; ui->pbProgress->setValue(val); //this->update(); } @@ -2084,9 +2085,9 @@ void FormToolbox::on_btnGeoColor_clicked() QApplication::setOverrideCursor(Qt::WaitCursor); // restore in FormImage::slotUpdateGeosat() ui->pbProgress->reset(); - if(geoindex == (int)eGeoSatellite::MET_11 || geoindex == (int)eGeoSatellite::MET_8) + if(geoindex == (int)eGeoSatellite::MET_11 || geoindex == (int)eGeoSatellite::MET_9 || geoindex == (int)eGeoSatellite::MET_8) ui->pbProgress->setMaximum(8+8+8); - else if(geoindex == (int)eGeoSatellite::MET_9 || geoindex == (int)eGeoSatellite::MET_10) + else if(geoindex == (int)eGeoSatellite::MET_10) ui->pbProgress->setMaximum(3+3+3); else if(geoindex == (int)eGeoSatellite::GOMS2) ui->pbProgress->setMaximum(6+6+6); @@ -2109,24 +2110,31 @@ void FormToolbox::on_btnRecipes_clicked() if(ui->lstRGB->currentRow() == -1) return; - if(!(geoindex == (int)eGeoSatellite::MET_11 || geoindex == (int)eGeoSatellite::MET_10 || geoindex == (int)eGeoSatellite::MET_8)) + if(!(geoindex == (int)eGeoSatellite::MET_11 || geoindex == (int)eGeoSatellite::MET_10 || geoindex == (int)eGeoSatellite::MET_9 || geoindex == (int)eGeoSatellite::MET_8)) return; QApplication::setOverrideCursor(Qt::WaitCursor); // restore in FormImage::slotUpdateGeosat() + formimage->displayImage(IMAGE_GEOSTATIONARY); + formimage->adjustPicSize(true); + ui->pbProgress->reset(); - ui->pbProgress->setMaximum(opts.geosatellites.at(geoindex).maxsegments); + ui->pbProgress->setMaximum(100); segs->seglgeo[geoindex]->areatype = 0; segs->seglgeo[geoindex]->setKindofImage("VIS_IR"); formimage->setKindOfImage("VIS_IR"); - formimage->displayImage(IMAGE_GEOSTATIONARY); - formimage->adjustPicSize(true); setToolboxButtons(false); - emit getrgbrecipe(ui->lstRGB->currentRow()); + emit switchstackedwidget(3); + + emit creatergbrecipe(ui->lstRGB->currentRow()); + + formimage->displayImage(IMAGE_GEOSTATIONARY); + formimage->adjustPicSize(true); + setToolboxButtons(true); QApplication::restoreOverrideCursor(); // restore in FormImage::slotUpdateGeosat() @@ -2184,7 +2192,7 @@ void FormToolbox::on_btnHRV_clicked() ui->pbProgress->reset(); - if(geoindex == (int)eGeoSatellite::MET_11 || geoindex == (int)eGeoSatellite::MET_10 || geoindex == (int)eGeoSatellite::MET_8) + if(geoindex == (int)eGeoSatellite::MET_11 || geoindex == (int)eGeoSatellite::MET_9 || geoindex == (int)eGeoSatellite::MET_8) { if(ui->cmbHRVtype->currentIndex() == 0 && ui->chkColorHRV->isChecked() == false) ui->pbProgress->setMaximum(5); @@ -2196,7 +2204,7 @@ void FormToolbox::on_btnHRV_clicked() ui->pbProgress->setMaximum(8+8+8+24); } - if(geoindex == (int)eGeoSatellite::MET_9) + if(geoindex == (int)eGeoSatellite::MET_10) { if(ui->cmbHRVtype->currentIndex() == 0 && ui->chkColorHRV->isChecked() == false) ui->pbProgress->setMaximum(5); @@ -2257,7 +2265,7 @@ void FormToolbox::onButtonColorHRV(QString type) formimage->setKindOfImage(type); formimage->displayImage(IMAGE_GEOSTATIONARY); - formimage->adjustPicSize(true); + formimage->adjustPicSize(false); setToolboxButtons(false); @@ -2611,8 +2619,6 @@ void FormToolbox::on_tabWidget_currentChanged(int index) { qDebug() << "on_tabWidget_currentChanged(int index) index = " << index; -// if(ui->tabWidgetVIIRS->currentIndex() == 0) -// qDebug() << "begin ui->tabWidgetVIIRS->currentIndex() == 0"; if(!forminfrascales->isHidden()) forminfrascales->hide(); diff --git a/core/formtoolbox.h b/core/formtoolbox.h index b4529b7..7785b35 100644 --- a/core/formtoolbox.h +++ b/core/formtoolbox.h @@ -58,6 +58,8 @@ class FormToolbox : public QWidget void setOLCIefrConfigsettings(); void setSLSTRConfigsettings(); void setComboGeo(int geoindex); + void setValuePrgBar(int val); + eProjectionType currentProjectionType; @@ -122,7 +124,8 @@ public slots: void getgeosatchannel(QString, QVector, QVector, int, bool); void overlaycorrection(int,int); void switchstackedwidget(int); - void getrgbrecipe(int recipe); + void creatergbrecipe(int recipe); + private slots: void on_btnCol_clicked(); diff --git a/core/formtoolbox.ui b/core/formtoolbox.ui index 1a9cff8..7ce6cab 100644 --- a/core/formtoolbox.ui +++ b/core/formtoolbox.ui @@ -5139,7 +5139,7 @@ - false + true QFrame::StyledPanel diff --git a/core/generalverticalperspective.cpp b/core/generalverticalperspective.cpp index 88fa6be..1604a95 100644 --- a/core/generalverticalperspective.cpp +++ b/core/generalverticalperspective.cpp @@ -295,8 +295,8 @@ void GeneralVerticalPerspective::CreateMapFromGeoStationary() { if(pixconv.geocoord2pixcoord(sub_lon, lat_rad*180.0/PI, lon_rad*180.0/PI, sl->COFF, sl->LOFF, sl->CFAC, sl->LFAC, &col, &row) == 0) { - row+=5; - col+=3; + row+=0; //5; + col+=50; //3; picrow = row; if( hrvmap == 0) { @@ -310,7 +310,7 @@ void GeneralVerticalPerspective::CreateMapFromGeoStationary() } else { - if(sl->getGeoSatellite() == eGeoSatellite::MET_9 || sl->getGeoSatellite() == eGeoSatellite::MET_10) + if(sl->getGeoSatellite() == eGeoSatellite::MET_10) { if( row < 5*464) { @@ -330,7 +330,7 @@ void GeneralVerticalPerspective::CreateMapFromGeoStationary() } } - else if(sl->getGeoSatellite() == eGeoSatellite::MET_11 || sl->getGeoSatellite() == eGeoSatellite::MET_8) + else if(sl->getGeoSatellite() == eGeoSatellite::MET_11 || sl->getGeoSatellite() == eGeoSatellite::MET_9 || sl->getGeoSatellite() == eGeoSatellite::MET_8) { if( row < (sl->areatype == 0 ? 5*464 : 11136)) { diff --git a/core/globals.h b/core/globals.h index 3bc32f8..763d1c6 100644 --- a/core/globals.h +++ b/core/globals.h @@ -80,6 +80,27 @@ enum class eGeoSatellite { NOGEO = 10 }; +//1 "Airmass RGB" +//2 "Dust RGB" +//3 "24 hours Microphysics RGB" +//4 "Ash RGB" +//5 "Day Microphysics RGB" +//6 "Severe Storms RGB" +//7 "Snow RGB" +//8 "Natural Colors RGB" +//9 "Night Microphysics RGB"; +enum class eRgbRecipes { + RGB_AIRMASS = 0, + RGB_DUST = 1, + RGB_24MICRO = 2, + RGB_ASH = 3, + RGB_DAYMICRO = 4, + RGB_STORMS = 5, + RGB_SNOW = 6, + RGB_NATURAL = 7, + RGB_NIGHTMICRO = 8 +}; + // structures struct floatVector { float x, y, z; diff --git a/core/globe.cpp b/core/globe.cpp index f00dc0f..fb9fbb6 100644 --- a/core/globe.cpp +++ b/core/globe.cpp @@ -1584,7 +1584,7 @@ void Globe::Render3DGeoSegmentLineFBO(int heightinimage, eGeoSatellite geo) rgbval = scanl[pix]; - if(geo == eGeoSatellite::MET_10 || geo == eGeoSatellite::MET_11) + if(geo == eGeoSatellite::MET_9 || geo == eGeoSatellite::MET_11) { if(pixconv.pixcoord2geocoord(segs->seglgeo[0]->geosatlon, pix, heightinimage, COFF_NONHRV, LOFF_NONHRV, CFAC_NONHRV, LFAC_NONHRV, &lat_deg, &lon_deg) == 0) { @@ -1597,7 +1597,7 @@ void Globe::Render3DGeoSegmentLineFBO(int heightinimage, eGeoSatellite geo) rainbow.append(qBlue(rgbval)); } } - else if(geo == eGeoSatellite::MET_9) + else if(geo == eGeoSatellite::MET_10) { if(pixconv.pixcoord2geocoord(segs->seglgeo[1]->geosatlon, pix, heightinimage, COFF_NONHRV, LOFF_NONHRV, CFAC_NONHRV, LFAC_NONHRV, &lat_deg, &lon_deg) == 0) { diff --git a/core/internal.cpp b/core/internal.cpp index 0f967b3..2c726ec 100644 --- a/core/internal.cpp +++ b/core/internal.cpp @@ -31,13 +31,13 @@ const std::string channel_spectrum[SEVIRI_N_BANDS] = const double band_solar_irradiance[][SEVIRI_N_BANDS] = { /* MSG-1, 321 */ - {65.2296, 73.0127, 62.3715, -999., -999., -999., -999., -999., -999., -999., 78.7599}, + {65.2296, 73.0127, 62.3715, -999., -999., -999., -999., -999., -999., -999., -999., 78.7599}, /* MSG-2, 322 */ - {65.2065, 73.1869, 61.9923, -999., -999., -999., -999., -999., -999., -999., 79.0113}, + {65.2065, 73.1869, 61.9923, -999., -999., -999., -999., -999., -999., -999., -999., 79.0113}, /* MSG-3, 323 */ - {65.5148, 73.1807, 62.0208, -999., -999., -999., -999., -999., -999., -999., 78.9416}, + {65.5148, 73.1807, 62.0208, -999., -999., -999., -999., -999., -999., -999., -999., 78.9416}, /* MSG-4, 324 */ - {65.2656, 73.1692, 61.9416, -999., -999., -999., -999., -999., -999., -999., 79.0035} + {65.2656, 73.1692, 61.9416, -999., -999., -999., -999., -999., -999., -999., -999., 79.0035} }; /* Ref: PDF_EFFECT_RAD_TO_BRIGHTNESS-1 */ diff --git a/core/lambertconformalconic.cpp b/core/lambertconformalconic.cpp index e90cc5a..365e661 100644 --- a/core/lambertconformalconic.cpp +++ b/core/lambertconformalconic.cpp @@ -386,7 +386,7 @@ void LambertConformalConic::CreateMapFromGeostationary() } else { - if(sl->getGeoSatellite() == eGeoSatellite::MET_9 || sl->getGeoSatellite() == eGeoSatellite::MET_10) + if(sl->getGeoSatellite() == eGeoSatellite::MET_10) { if( picrow >= 0 && picrow < 5*464) { diff --git a/core/main.cpp b/core/main.cpp index 529a2ce..d06ae8e 100644 --- a/core/main.cpp +++ b/core/main.cpp @@ -12,7 +12,7 @@ #include -#define APPVERSION "1.3.7" +#define APPVERSION "1.3.9" using namespace std; @@ -32,7 +32,6 @@ bool ptrimagebusy; void myMessageOutput(QtMsgType type, const QMessageLogContext &context, const QString &msg) { - #ifdef NDEBUG // release mode code QString strout; diff --git a/core/mainwindow.cpp b/core/mainwindow.cpp index 4ab89ac..a3318e9 100644 --- a/core/mainwindow.cpp +++ b/core/mainwindow.cpp @@ -85,6 +85,7 @@ MainWindow::MainWindow(QWidget *parent) : connect(seglist->seglM01hrpt, SIGNAL(progressCounter(int)), formtoolbox, SLOT(setValueProgressBar(int))); connect(seglist->seglM02hrpt, SIGNAL(progressCounter(int)), formtoolbox, SLOT(setValueProgressBar(int))); + connect(seglist->seglviirsdnb, SIGNAL(displayDNBGraph()), formtoolbox, SLOT(slotDisplayDNBGraph())); connect(seglist->seglviirsdnbnoaa20, SIGNAL(displayDNBGraph()), formtoolbox, SLOT(slotDisplayDNBGraph())); @@ -171,7 +172,7 @@ MainWindow::MainWindow(QWidget *parent) : connect( formtoolbox, SIGNAL(getgeosatchannel(QString, QVector, QVector, int, bool)), formgeostationary, SLOT(slotCreateGeoImage(QString, QVector, QVector, int, bool))); connect( formtoolbox, SIGNAL(switchstackedwidget(int)), this, SLOT(slotSwitchStackedWindow(int))); - connect( formtoolbox, SIGNAL(getrgbrecipe(int)), formgeostationary, SLOT(slotCreateRGBrecipe(int))); + connect( formtoolbox, SIGNAL(creatergbrecipe(int)), formgeostationary, SLOT(slotCreateRGBrecipe(int))); connect( formgeostationary, SIGNAL(enabletoolboxbuttons(bool)), formtoolbox, SLOT(setToolboxButtons(bool))); diff --git a/core/options.cpp b/core/options.cpp index 6f53b4a..4bed1eb 100644 --- a/core/options.cpp +++ b/core/options.cpp @@ -787,8 +787,8 @@ void Options::CreateGeoSatelliteIni() geosatellites[2].imageheight = 1392; geosatellites[2].imagewidthhrv0 = 5568; geosatellites[2].imageheighthrv0 = 2320; - geosatellites[2].imagewidthhrv1 = 0; - geosatellites[2].imageheighthrv1 = 0; + geosatellites[2].imagewidthhrv1 = 5568; + geosatellites[2].imageheighthrv1 = 11136; geosatellites[2].indexspectrum = 26; geosatellites[2].indexfilenbr = 36; @@ -811,8 +811,8 @@ void Options::CreateGeoSatelliteIni() geosatellites[2].segmentlengthhrv = 464; geosatellites[2].startsegmentnbrtype0 = 5; geosatellites[2].startsegmentnbrhrvtype0 = 19; - geosatellites[2].startsegmentnbrtype1 = 0; - geosatellites[2].startsegmentnbrhrvtype1 = 0; + geosatellites[2].startsegmentnbrtype1 = 1; + geosatellites[2].startsegmentnbrhrvtype1 = 1; geosatellites[2].prologfile = true; geosatellites[2].epilogfile = true; diff --git a/core/segmentimage.cpp b/core/segmentimage.cpp index 40d5470..b454d80 100644 --- a/core/segmentimage.cpp +++ b/core/segmentimage.cpp @@ -57,62 +57,89 @@ SegmentImage::SegmentImage() SetupRGBrecipes(); } +//void SegmentImage::SetupRGBrecipes() +//{ +// QList recipes; + +// recipes << "Airmass RGB" +// << "Dust RGB" +// << "24 hours Microphysics RGB" +// << "Ash RGB" +// << "Day Microphysics RGB" +// << "Severe Storms RGB" +// << "Snow RGB" +// << "Natural Colors RGB" +// << "Night Microphysics RGB"; + +//} void SegmentImage::SetupRGBrecipes() { + //0 "Airmass RGB" + //1 "Dust RGB" + //2 "24 hours Microphysics RGB" + //3 "Ash RGB" + //4 "Day Microphysics RGB" + //5 "Severe Storms RGB" + //6 "Snow RGB" + //7 "Natural Colors RGB" + //8 "Night Microphysics RGB"; + //*********** //Airmass RGB //*********** RGBRecipe airmass; airmass.Name = "Airmass RGB"; - RGBRecipeColor Red; - RGBRecipeColor Green; - RGBRecipeColor Blue; - - Red.channels.append("WV_062"); - Red.channels.append("WV_073"); - Red.spectral_channel_nbr.append(GetSpectralChannelNbr("WV_062")); - Red.spectral_channel_nbr.append(GetSpectralChannelNbr("WV_073")); - Red.subtract.append(false); - Red.subtract.append(true); - Red.inverse.append(false); - Red.inverse.append(false); - Red.reflective.append(false); - Red.reflective.append(false); - Red.rangefrom = -25.0; - Red.rangeto = 0.0; - Red.dimension = "K"; - Red.gamma = 1.0; - Red.units = SEVIRI_UNIT_BT; - airmass.Colorvector.append(Red); - - Green.channels.append("IR_097"); - Green.channels.append("IR_108"); - Green.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_097")); - Green.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); - Green.subtract.append(false); - Green.subtract.append(true); - Green.inverse.append(false); - Green.inverse.append(false); - Green.reflective.append(false); - Green.reflective.append(false); - Green.rangefrom = -40.0; - Green.rangeto = 5.0; - Green.dimension = "K"; - Green.gamma = 1.0; - Green.units = SEVIRI_UNIT_BT; - airmass.Colorvector.append(Green); - - Blue.channels.append("WV_062"); - Blue.spectral_channel_nbr.append(GetSpectralChannelNbr("WV_062")); - Blue.subtract.append(false); - Blue.inverse.append(true); - Blue.reflective.append(false); - Blue.rangefrom = 243.0; - Blue.rangeto = 208.0; - Blue.dimension = "K"; - Blue.gamma = 1.0; - Blue.units = SEVIRI_UNIT_BT; - airmass.Colorvector.append(Blue); + airmass.reflectivechannel = false; + + RGBRecipeColor Red1; + RGBRecipeColor Green1; + RGBRecipeColor Blue1; + + Red1.channels.append("WV_062"); + Red1.channels.append("WV_073"); + Red1.spectral_channel_nbr.append(GetSpectralChannelNbr("WV_062")); + Red1.spectral_channel_nbr.append(GetSpectralChannelNbr("WV_073")); + Red1.subtract.append(false); + Red1.subtract.append(true); + Red1.inverse.append(false); + Red1.inverse.append(false); + Red1.reflective.append(false); + Red1.reflective.append(false); + Red1.rangefrom = -25.0; + Red1.rangeto = 0.0; + Red1.dimension = "K"; + Red1.gamma = 1.0; + Red1.units = SEVIRI_UNIT_BT; + airmass.Colorvector.append(Red1); + + Green1.channels.append("IR_097"); + Green1.channels.append("IR_108"); + Green1.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_097")); + Green1.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); + Green1.subtract.append(false); + Green1.subtract.append(true); + Green1.inverse.append(false); + Green1.inverse.append(false); + Green1.reflective.append(false); + Green1.reflective.append(false); + Green1.rangefrom = -40.0; + Green1.rangeto = 5.0; + Green1.dimension = "K"; + Green1.gamma = 1.0; + Green1.units = SEVIRI_UNIT_BT; + airmass.Colorvector.append(Green1); + + Blue1.channels.append("WV_062"); + Blue1.spectral_channel_nbr.append(GetSpectralChannelNbr("WV_062")); + Blue1.subtract.append(false); + Blue1.inverse.append(true); + Blue1.reflective.append(false); + Blue1.rangefrom = 208.0; + Blue1.rangeto = 243.0; + Blue1.dimension = "K"; + Blue1.gamma = 1.0; + Blue1.units = SEVIRI_UNIT_BT; + airmass.Colorvector.append(Blue1); rgbrecipes.append(airmass); @@ -121,6 +148,8 @@ void SegmentImage::SetupRGBrecipes() //********* RGBRecipe dust; dust.Name = "Dust RGB"; + dust.reflectivechannel = false; + RGBRecipeColor Red2; RGBRecipeColor Green2; RGBRecipeColor Blue2; @@ -142,7 +171,6 @@ void SegmentImage::SetupRGBrecipes() Red2.units = SEVIRI_UNIT_BT; dust.Colorvector.append(Red2); - Green2.channels.clear(); Green2.channels.append("IR_108"); Green2.channels.append("IR_087"); Green2.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); @@ -156,11 +184,10 @@ void SegmentImage::SetupRGBrecipes() Green2.rangefrom = 0.0; Green2.rangeto = 15.0; Green2.dimension = "K"; - Green2.gamma = 5.0; //2.5; + Green2.gamma = 2.5; Green2.units = SEVIRI_UNIT_BT; dust.Colorvector.append(Green2); - Blue2.channels.clear(); Blue2.channels.append("IR_108"); Blue2.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); Blue2.subtract.append(false); @@ -176,116 +203,460 @@ void SegmentImage::SetupRGBrecipes() rgbrecipes.append(dust); //****************** - //Natual Colours RGB + //24 hours Microphysics RGB //****************** - RGBRecipe natural; - natural.Name = "Natural Colours RGB"; + + RGBRecipe micro24; + micro24.Name = "24 hours Microphysics RGB"; + micro24.reflectivechannel = false; + RGBRecipeColor Red3; RGBRecipeColor Green3; RGBRecipeColor Blue3; - Red3.channels.append("IR_016"); - Red3.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_016")); + Red3.channels.append("IR_120"); + Red3.channels.append("IR_108"); + Red3.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_120")); + Red3.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); Red3.subtract.append(false); + Red3.subtract.append(true); + Red3.inverse.append(false); Red3.inverse.append(false); Red3.reflective.append(false); - Red3.rangefrom = 0.0; - Red3.rangeto = 100.0; - Red3.dimension = "%"; + Red3.reflective.append(false); + Red3.rangefrom = -4.0; + Red3.rangeto = +2.0; + Red3.dimension = "K"; Red3.gamma = 1.0; - Red3.units = SEVIRI_UNIT_REF; - natural.Colorvector.append(Red3); + Red3.units = SEVIRI_UNIT_BT; + micro24.Colorvector.append(Red3); - Green3.channels.clear(); - Green3.channels.append("VIS008"); - Green3.spectral_channel_nbr.append(GetSpectralChannelNbr("VIS008")); + Green3.channels.append("IR_108"); + Green3.channels.append("IR_087"); + Green3.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); + Green3.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_087")); Green3.subtract.append(false); + Green3.subtract.append(true); Green3.inverse.append(false); + Green3.inverse.append(false); + Green3.reflective.append(false); Green3.reflective.append(false); Green3.rangefrom = 0.0; - Green3.rangeto = 100.0; - Green3.dimension = "%"; - Green3.gamma = 1.0; - Green3.units = SEVIRI_UNIT_REF; - natural.Colorvector.append(Green3); - - Blue3.channels.clear(); - Blue3.channels.append("VIS006"); - Blue3.spectral_channel_nbr.append(GetSpectralChannelNbr("VIS006")); + Green3.rangeto = 6.0; + Green3.dimension = "K"; + Green3.gamma = 1.2; + Green3.units = SEVIRI_UNIT_BT; + micro24.Colorvector.append(Green3); + + Blue3.channels.append("IR_108"); + Blue3.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); Blue3.subtract.append(false); Blue3.inverse.append(false); Blue3.reflective.append(false); - Blue3.rangefrom = 0.0; - Blue3.rangeto = 100.0; - Blue3.dimension = "%"; + Blue3.rangefrom = 248.0; + Blue3.rangeto = 303.0; + Blue3.dimension = "K"; Blue3.gamma = 1.0; - Blue3.units = SEVIRI_UNIT_REF; - natural.Colorvector.append(Blue3); + Blue3.units = SEVIRI_UNIT_BT; + micro24.Colorvector.append(Blue3); - rgbrecipes.append(natural); + rgbrecipes.append(micro24); //****************** - //Test + // Ash RGB //****************** - RGBRecipe test; - test.Name = "Test RGB"; + + RGBRecipe ash; + ash.Name = "Ash RGB"; + ash.reflectivechannel = false; + RGBRecipeColor Red4; RGBRecipeColor Green4; RGBRecipeColor Blue4; - Red4.channels.append("IR_087"); - Red4.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_087")); + Red4.channels.append("IR_120"); + Red4.channels.append("IR_108"); + Red4.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_120")); + Red4.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); Red4.subtract.append(false); + Red4.subtract.append(true); + Red4.inverse.append(false); Red4.inverse.append(false); Red4.reflective.append(false); - Red4.rangefrom = 0.0; - Red4.rangeto = 400.0; + Red4.reflective.append(false); + Red4.rangefrom = -4.0; + Red4.rangeto = +2.0; Red4.dimension = "K"; Red4.gamma = 1.0; Red4.units = SEVIRI_UNIT_BT; - test.Colorvector.append(Red4); + ash.Colorvector.append(Red4); - Green4.channels.clear(); + Green4.channels.append("IR_108"); Green4.channels.append("IR_087"); + Green4.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); Green4.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_087")); Green4.subtract.append(false); + Green4.subtract.append(true); + Green4.inverse.append(false); Green4.inverse.append(false); Green4.reflective.append(false); - Green4.rangefrom = 0.0; - Green4.rangeto = 400.0; + Green4.reflective.append(false); + Green4.rangefrom = -4.0; + Green4.rangeto = 5.0; Green4.dimension = "K"; Green4.gamma = 1.0; Green4.units = SEVIRI_UNIT_BT; - test.Colorvector.append(Green4); + ash.Colorvector.append(Green4); - Blue4.channels.clear(); - Blue4.channels.append("IR_087"); - Blue4.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_087")); + Blue4.channels.append("IR_108"); + Blue4.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); Blue4.subtract.append(false); Blue4.inverse.append(false); Blue4.reflective.append(false); - Blue4.rangefrom = 0.0; - Blue4.rangeto = 400.0; + Blue4.rangefrom = 243.0; + Blue4.rangeto = 303.0; Blue4.dimension = "K"; Blue4.gamma = 1.0; Blue4.units = SEVIRI_UNIT_BT; - test.Colorvector.append(Blue4); + ash.Colorvector.append(Blue4); - rgbrecipes.append(test); + rgbrecipes.append(ash); - int recipe = 0; - for(int colorindex = 0; colorindex < 3; colorindex++) + //****************** + //Day Microphysics RGB + //****************** + RGBRecipe day; + day.Name = "Day Microphysics RGB"; + day.reflectivechannel = true; + + RGBRecipeColor Red5; + RGBRecipeColor Green5; + RGBRecipeColor Blue5; + + Red5.channels.append("VIS008"); + Red5.spectral_channel_nbr.append(GetSpectralChannelNbr("VIS008")); + Red5.subtract.append(false); + Red5.inverse.append(false); + Red5.reflective.append(false); + Red5.rangefrom = 0.0; + Red5.rangeto = 1.0; + Red5.dimension = "%"; + Red5.gamma = 1.0; + Red5.units = SEVIRI_UNIT_BRF; + day.Colorvector.append(Red5); + + Green5.channels.append("IR_039"); + Green5.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_039")); + Green5.subtract.append(false); + Green5.inverse.append(false); + Green5.reflective.append(true); + Green5.rangefrom = 0.0; + Green5.rangeto = 1.0; + Green5.dimension = "%"; + Green5.gamma = 2.5; + Green5.units = SEVIRI_UNIT_REFL39; + day.Colorvector.append(Green5); + + Blue5.channels.append("IR_108"); + Blue5.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); + Blue5.subtract.append(false); + Blue5.inverse.append(false); + Blue5.reflective.append(false); + Blue5.rangefrom = 203.0; + Blue5.rangeto = 323.0; + Blue5.dimension = "K"; + Blue5.gamma = 1.0; + Blue5.units = SEVIRI_UNIT_BT; + day.Colorvector.append(Blue5); + + rgbrecipes.append(day); + + //****************** + // Severe Storms RGB + //****************** + + RGBRecipe storm; + storm.Name = "Severe Storms RGB"; + storm.reflectivechannel = false; + + RGBRecipeColor Red6; + RGBRecipeColor Green6; + RGBRecipeColor Blue6; + + Red6.channels.append("WV_062"); + Red6.channels.append("WV_073"); + Red6.spectral_channel_nbr.append(GetSpectralChannelNbr("WV_062")); + Red6.spectral_channel_nbr.append(GetSpectralChannelNbr("WV_073")); + Red6.subtract.append(false); + Red6.subtract.append(true); + Red6.inverse.append(false); + Red6.inverse.append(false); + Red6.reflective.append(false); + Red6.reflective.append(false); + Red6.rangefrom = -35.0; + Red6.rangeto = +5.0; + Red6.dimension = "K"; + Red6.gamma = 1.0; + Red6.units = SEVIRI_UNIT_BT; + storm.Colorvector.append(Red6); + + Green6.channels.append("IR_039"); + Green6.channels.append("IR_108"); + Green6.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_039")); + Green6.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); + Green6.subtract.append(false); + Green6.subtract.append(true); + Green6.inverse.append(false); + Green6.inverse.append(false); + Green6.reflective.append(false); + Green6.reflective.append(false); + Green6.rangefrom = -5.0; + Green6.rangeto = 60.0; + Green6.dimension = "K"; + Green6.gamma = 0.5; + Green6.units = SEVIRI_UNIT_BT; + storm.Colorvector.append(Green6); + + Blue6.channels.append("IR_016"); + Blue6.channels.append("VIS006"); + Blue6.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_016")); + Blue6.spectral_channel_nbr.append(GetSpectralChannelNbr("VIS006")); + Blue6.subtract.append(false); + Blue6.subtract.append(true); + Blue6.inverse.append(false); + Blue6.inverse.append(false); + Blue6.reflective.append(false); + Blue6.reflective.append(false); + Blue6.rangefrom = -75.0; + Blue6.rangeto = 25.0; + Blue6.dimension = "%"; + Blue6.gamma = 1.0; + Blue6.units = SEVIRI_UNIT_REF; + storm.Colorvector.append(Blue6); + + rgbrecipes.append(storm); + + //****************** + //Snow RGB + //****************** + RGBRecipe snow; + snow.Name = "Snow RGB"; + snow.reflectivechannel = true; + + RGBRecipeColor Red7; + RGBRecipeColor Green7; + RGBRecipeColor Blue7; + + Red7.channels.append("VIS008"); + Red7.spectral_channel_nbr.append(GetSpectralChannelNbr("VIS008")); + Red7.subtract.append(false); + Red7.inverse.append(false); + Red7.reflective.append(false); + Red7.rangefrom = 0.0; + Red7.rangeto = 1.0; + Red7.dimension = "%"; + Red7.gamma = 1.0; + Red7.units = SEVIRI_UNIT_REF; + snow.Colorvector.append(Red7); + + Green7.channels.append("IR_016"); + Green7.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_016")); + Green7.subtract.append(false); + Green7.inverse.append(false); + Green7.reflective.append(false); + Green7.rangefrom = 0.0; + Green7.rangeto = 0.7; + Green7.dimension = "%"; + Green7.gamma = 1.7; + Green7.units = SEVIRI_UNIT_REF; + snow.Colorvector.append(Green7); + + Blue7.channels.append("IR_039"); + Blue7.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_039")); + Blue7.subtract.append(false); + Blue7.inverse.append(false); + Blue7.reflective.append(true); + Blue7.rangefrom = 0.0; + Blue7.rangeto = 0.3; + Blue7.dimension = "%"; + Blue7.gamma = 1.7; + Blue7.units = SEVIRI_UNIT_REF; + snow.Colorvector.append(Blue7); + + rgbrecipes.append(snow); + + //****************** + //Natual Colours RGB + //****************** + RGBRecipe natural; + natural.Name = "Natural Colours RGB"; + natural.reflectivechannel = false; + + RGBRecipeColor Red8; + RGBRecipeColor Green8; + RGBRecipeColor Blue8; + + Red8.channels.append("IR_016"); + Red8.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_016")); + Red8.subtract.append(false); + Red8.inverse.append(false); + Red8.reflective.append(false); + Red8.rangefrom = 0.0; + Red8.rangeto = 1.0; + Red8.dimension = "%"; + Red8.gamma = 1.0; + Red8.units = SEVIRI_UNIT_REF; + natural.Colorvector.append(Red8); + + Green8.channels.append("VIS008"); + Green8.spectral_channel_nbr.append(GetSpectralChannelNbr("VIS008")); + Green8.subtract.append(false); + Green8.inverse.append(false); + Green8.reflective.append(false); + Green8.rangefrom = 0.0; + Green8.rangeto = 1.0; + Green8.dimension = "%"; + Green8.gamma = 1.0; + Green8.units = SEVIRI_UNIT_REF; + natural.Colorvector.append(Green8); + + Blue8.channels.append("VIS006"); + Blue8.spectral_channel_nbr.append(GetSpectralChannelNbr("VIS006")); + Blue8.subtract.append(false); + Blue8.inverse.append(false); + Blue8.reflective.append(false); + Blue8.rangefrom = 0.0; + Blue8.rangeto = 1.0; + Blue8.dimension = "%"; + Blue8.gamma = 1.0; + Blue8.units = SEVIRI_UNIT_REF; + natural.Colorvector.append(Blue8); + + rgbrecipes.append(natural); + + //*********** + //Night Microphysics RGB + //*********** + RGBRecipe night; + night.Name = "Night Microphysics RGB"; + night.reflectivechannel = false; + + RGBRecipeColor Red9; + RGBRecipeColor Green9; + RGBRecipeColor Blue9; + + Red9.channels.append("IR_120"); + Red9.channels.append("IR_108"); + Red9.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_120")); + Red9.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); + Red9.subtract.append(false); + Red9.subtract.append(true); + Red9.inverse.append(false); + Red9.inverse.append(false); + Red9.reflective.append(false); + Red9.reflective.append(false); + Red9.rangefrom = -4.0; + Red9.rangeto = 2.0; + Red9.dimension = "K"; + Red9.gamma = 1.0; + Red9.units = SEVIRI_UNIT_BT; + night.Colorvector.append(Red9); + + Green9.channels.append("IR_108"); + Green9.channels.append("IR_039"); + Green9.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); + Green9.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_039")); + Green9.subtract.append(false); + Green9.subtract.append(true); + Green9.inverse.append(false); + Green9.inverse.append(false); + Green9.reflective.append(false); + Green9.reflective.append(false); + Green9.rangefrom = 0.0; + Green9.rangeto = 10.0; + Green9.dimension = "K"; + Green9.gamma = 1.0; + Green9.units = SEVIRI_UNIT_BT; + night.Colorvector.append(Green9); + + Blue9.channels.append("IR_108"); + Blue9.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_108")); + Blue9.subtract.append(false); + Blue9.inverse.append(false); + Blue9.reflective.append(false); + Blue9.rangefrom = 243.0; + Blue9.rangeto = 293.0; + Blue9.dimension = "K"; + Blue9.gamma = 1.0; + Blue9.units = SEVIRI_UNIT_BT; + night.Colorvector.append(Blue9); + + rgbrecipes.append(night); + + + //****************** + //Test RGB + //****************** +// RGBRecipe test; +// test.Name = "Test RGB"; +// test.reflectivechannel = false; + +// RGBRecipeColor Red10; +// RGBRecipeColor Green10; +// RGBRecipeColor Blue10; + +// Red10.channels.append("IR_039"); +// Red10.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_039")); +// Red10.subtract.append(false); +// Red10.inverse.append(false); +// Red10.reflective.append(false); +// Red10.rangefrom = 250.0; +// Red10.rangeto = 355.55; +// Red10.dimension = "K"; +// Red10.gamma = 1.0; +// Red10.units = SEVIRI_UNIT_BT; +// test.Colorvector.append(Red10); + +// Green10.channels.append("IR_039"); +// Green10.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_039")); +// Green10.subtract.append(false); +// Green10.inverse.append(false); +// Green10.reflective.append(false); +// Green10.rangefrom = 250.0; +// Green10.rangeto = 355.55; +// Green10.dimension = "K"; +// Green10.gamma = 1.0; +// Green10.units = SEVIRI_UNIT_BT; +// test.Colorvector.append(Green10); + +// Blue10.channels.append("IR_039"); +// Blue10.spectral_channel_nbr.append(GetSpectralChannelNbr("IR_039")); +// Blue10.subtract.append(false); +// Blue10.inverse.append(false); +// Blue10.reflective.append(false); +// Blue10.rangefrom = 250.0; +// Blue10.rangeto = 355.55; +// Blue10.dimension = "K"; +// Blue10.gamma = 1.0; +// Blue10.units = SEVIRI_UNIT_BT; +// test.Colorvector.append(Blue10); + +// rgbrecipes.append(test); + + for( int recipe = 0; recipe < rgbrecipes.size(); recipe++) { - for(int i = 0; i < rgbrecipes[recipe].Colorvector.at(colorindex).channels.length(); i++) + for(int colorindex = 0; colorindex < 3; colorindex++) { - qDebug() << colorindex << " " << rgbrecipes[recipe].Colorvector.at(colorindex).channels.at(i) << " " << - rgbrecipes[recipe].Colorvector.at(colorindex).spectral_channel_nbr.at(i) << " " << - rgbrecipes[recipe].Colorvector.at(colorindex).subtract.at(i); + for(int i = 0; i < rgbrecipes[recipe].Colorvector.at(colorindex).channels.length(); i++) + { + qDebug() << rgbrecipes[recipe].Name << " " << colorindex << " " << rgbrecipes[recipe].Colorvector.at(colorindex).channels.at(i) << " " << + rgbrecipes[recipe].Colorvector.at(colorindex).spectral_channel_nbr.at(i) << " " << + rgbrecipes[recipe].Colorvector.at(colorindex).subtract.at(i); + } } } - } int SegmentImage::GetSpectralChannelNbr(QString channel) diff --git a/core/segmentimage.h b/core/segmentimage.h index 72d29ba..e7a2bfb 100644 --- a/core/segmentimage.h +++ b/core/segmentimage.h @@ -18,10 +18,27 @@ enum seviri_units { SEVIRI_UNIT_REF, SEVIRI_UNIT_BRF, SEVIRI_UNIT_BT, + SEVIRI_UNIT_REFL39, N_SEVIRI_UNITS }; +struct Bandstorage { + int listindex; + int spectral_channel_nbr; + QString directory; + QString productid1; + QString productid2; + QString timing; + double day_of_year; + float min; + float max; + float *data; + seviri_units units; + double slope; + double offset; +}; + struct RGBRecipeColor { QStringList channels; QList spectral_channel_nbr; // 1 - 12 @@ -38,6 +55,7 @@ struct RGBRecipeColor { struct RGBRecipe { QString Name; + bool reflectivechannel; QVector Colorvector; }; diff --git a/core/segmentlistgeostationary.cpp b/core/segmentlistgeostationary.cpp index 2ca74e1..27f319b 100644 --- a/core/segmentlistgeostationary.cpp +++ b/core/segmentlistgeostationary.cpp @@ -2,9 +2,10 @@ #include #include #include "segmentlistgeostationary.h" -#include "segmentimage.h" #include "qcompressor.h" #include "pixgeoconversion.h" +#include "misc_util.h" +#include "internal.h" #ifdef _WIN32 #include @@ -99,6 +100,10 @@ void SegmentListGeostationary::doComposeGeostationarynetCDFInThread(SegmentListG sm->ComposeSegmentImagenetCDFInThread(); //strlist, spectrumvector, inversevector, histogrammethod); } +void SegmentListGeostationary::doComposeGeoRGBRecipe(SegmentListGeostationary *sm, int recipe) +{ + sm->ComposeGeoRGBRecipeInThread(recipe); +} SegmentListGeostationary::SegmentListGeostationary(QObject *parent, int geoindex) : QObject(parent) @@ -589,7 +594,7 @@ void SegmentListGeostationary::ComposeSegmentImageXRIT( QString filepath, int ch emit signalcomposefinished(kindofimage); } - if((kindofimage == "VIS_IR" || kindofimage == "VIS_IR Color") && allSegmentsReceived()) + if((kindofimage == "VIS_IR" || kindofimage == "VIS_IR Color") && allVIS_IRSegmentsReceived()) { qDebug() << "-----> VIS_IR Color | VIS_IR | HRV and allSegmentsReceived"; this->ComposeVISIR(); @@ -742,7 +747,7 @@ void SegmentListGeostationary::ComposeSegmentImageXRITHimawari( QString filepath - if(allSegmentsReceived()) + if(allVIS_IRSegmentsReceived()) { ComposeVISIRHimawari(); emit signalcomposefinished(kindofimage); @@ -1766,6 +1771,8 @@ void SegmentListGeostationary::ComposeSegmentImagenetCDFInThread() //(QStringLis emit signalcomposefinished(kindofimage); } + + void SegmentListGeostationary::CalculateLUTGeo(int colorindex) { qDebug() << "start SegmentListGeostationary::CalculateLUTGeo()"; @@ -2454,10 +2461,10 @@ bool SegmentListGeostationary::allHRVColorSegmentsReceived() return true; } -bool SegmentListGeostationary::allSegmentsReceived() +bool SegmentListGeostationary::allVIS_IRSegmentsReceived() { - qDebug() << QString("SegmentListGeostationary::allSegmentsReceived()"); + qDebug() << QString("SegmentListGeostationary::allVIS_IRSegmentsReceived()"); int pbCounter = 0; int pbCounterRed = 0; @@ -2473,7 +2480,7 @@ bool SegmentListGeostationary::allSegmentsReceived() } else if(m_GeoSatellite == eGeoSatellite::MET_11 || m_GeoSatellite == eGeoSatellite::MET_10 || m_GeoSatellite == eGeoSatellite::MET_9 || m_GeoSatellite == eGeoSatellite::MET_8) { - for(int i = (m_GeoSatellite == eGeoSatellite::MET_9 ? 5 : 0) ; i < 8; i++) + for(int i = (this->bisRSS ? 5 : 0) ; i < 8; i++) { if (isPresentRed[i] && issegmentcomposedRed[i] == true) { @@ -2537,7 +2544,7 @@ bool SegmentListGeostationary::allSegmentsReceived() } else if(m_GeoSatellite == eGeoSatellite::MET_9) { - for(int i = 5; i < 8; i++) + for(int i = 0; i < 8; i++) { if (isPresentRed[i] && issegmentcomposedRed[i] == true) pbCounter++; @@ -2571,7 +2578,7 @@ bool SegmentListGeostationary::allSegmentsReceived() { if(m_GeoSatellite == eGeoSatellite::MET_11 || m_GeoSatellite == eGeoSatellite::MET_10 || m_GeoSatellite == eGeoSatellite::MET_9 || m_GeoSatellite == eGeoSatellite::MET_8) { - for(int i = (m_GeoSatellite == eGeoSatellite::MET_9 ? 5 : 0) ; i < 8; i++) + for(int i = (this->bisRSS ? 5 : 0) ; i < 8; i++) { if (isPresentRed[i] && issegmentcomposedRed[i] == true) pbCounter++; @@ -2593,7 +2600,7 @@ bool SegmentListGeostationary::allSegmentsReceived() } } - qDebug() << QString("SegmentListGeostationary::allSegmentsReceived() pbCounter = %1").arg(pbCounter); + qDebug() << QString("SegmentListGeostationary::allVIS_IRSegmentsReceived() pbCounter = %1").arg(pbCounter); if(pbCounterRed == 8) qDebug() << QString("pbCounterRed = 8"); @@ -2609,7 +2616,7 @@ bool SegmentListGeostationary::allSegmentsReceived() { if(m_GeoSatellite == eGeoSatellite::MET_11 || m_GeoSatellite == eGeoSatellite::MET_10 || m_GeoSatellite == eGeoSatellite::MET_9 || m_GeoSatellite == eGeoSatellite::MET_8) { - for(int i = (m_GeoSatellite == eGeoSatellite::MET_9 ? 5 : 0) ; i < 8; i++) + for(int i = (this->bisRSS ? 5 : 0) ; i < 8; i++) { if (isPresentRed[i] && issegmentcomposedRed[i] == false) return false; @@ -2655,7 +2662,7 @@ bool SegmentListGeostationary::allSegmentsReceived() } else if(m_GeoSatellite == eGeoSatellite::MET_9) { - for(int i = 5; i < 8; i++) + for(int i = 0; i < 8; i++) { if (isPresentRed[i] && issegmentcomposedRed[i] == false) return false; @@ -2716,7 +2723,7 @@ bool SegmentListGeostationary::allSegmentsReceived() } } - qDebug() << "SegmentListGeostationary::allSegmentsReceived() returns true"; + qDebug() << "SegmentListGeostationary::allVIS_IRSegmentsReceived() returns true"; emit progressCounter(100); return true; @@ -2882,5 +2889,1449 @@ void SegmentListGeostationary::SetupContrastStretch(quint16 x1, quint16 y1, quin //qDebug() << QString("A1 = %1;B1 = %2;A2 = %3;B2 = %4;A3 = %5;B3 = %6").arg(A1).arg(B1).arg(A2).arg(B2).arg(A3).arg(B3); } +void SegmentListGeostationary::ComposeGeoRGBRecipe(int recipe, QString tex) +{ + this->tex = tex; + QtConcurrent::run(doComposeGeoRGBRecipe, this, recipe); +} + +void SegmentListGeostationary::ComposeGeoRGBRecipeInThread(int recipe) +{ + qDebug() << "start SegmentListGeostationary::ComposeGeoRGBRecipeInThread(int recipe)"; + emit progressCounter(10); + + QStringList redbandlist = imageptrs->rgbrecipes[recipe].Colorvector.at(0).channels; + QStringList greenbandlist = imageptrs->rgbrecipes[recipe].Colorvector.at(1).channels; + QStringList bluebandlist = imageptrs->rgbrecipes[recipe].Colorvector.at(2).channels; + QList uniqueunitlist; + seviri_units redunits = imageptrs->rgbrecipes[recipe].Colorvector.at(0).units; + seviri_units greenunits = imageptrs->rgbrecipes[recipe].Colorvector.at(1).units; + seviri_units blueunits = imageptrs->rgbrecipes[recipe].Colorvector.at(2).units; + QList uniquechannelnbrlist; + + QStringList uniquebandlist; + + for(int i = 0; i < redbandlist.length(); i++) + { + if(!uniquebandlist.contains(redbandlist.at(i))) + { + uniquebandlist.append(redbandlist.at(i)); + uniqueunitlist.append(redunits); + uniquechannelnbrlist.append(imageptrs->rgbrecipes[recipe].Colorvector.at(0).spectral_channel_nbr.at(i)); + } + } + for(int i = 0; i < greenbandlist.length(); i++) + { + if(!uniquebandlist.contains(greenbandlist.at(i))) + { + uniquebandlist.append(greenbandlist.at(i)); + uniqueunitlist.append(greenunits); + uniquechannelnbrlist.append(imageptrs->rgbrecipes[recipe].Colorvector.at(1).spectral_channel_nbr.at(i)); + } + } + for(int i = 0; i < bluebandlist.length(); i++) + { + if(!uniquebandlist.contains(bluebandlist.at(i))) + { + uniquebandlist.append(bluebandlist.at(i)); + uniqueunitlist.append(blueunits); + uniquechannelnbrlist.append(imageptrs->rgbrecipes[recipe].Colorvector.at(2).spectral_channel_nbr.at(i)); + } + } + + + QString recipename = imageptrs->rgbrecipes[recipe].Name; + qDebug() << "recipename = " << recipename; + qDebug() << "uniquebandlist contains = " << uniquebandlist.count(); + for(int i = 0; i < uniquebandlist.length(); i++) + { + qDebug() << "channels at " << uniquechannelnbrlist.at(i) << " = " << uniquebandlist.at(i) << " unit = " << uniqueunitlist.at(i); + } + + + qDebug() << "tex = " << this->tex << " imagepath = " << this->getImagePath(); + //tex = "2018-02-13 14:45" + + this->ResetSegments(); + imageptrs->ResetPtrImage(); + + QString directory = this->getImagePath(); + QString productid1 = opts.geosatellites.at(geoindex).searchstring.mid(6, 4); + QString productid2 = "IR_016"; + QString timing = tex.mid(0, 4) + tex.mid(5, 2) + tex.mid(8, 2) + tex.mid(13, 2) + tex.mid(16, 2); + + qDebug() << "directory = " << directory; + qDebug() << "productid1 = " << productid1; + qDebug() << "productid2 = " << productid2; + qDebug() << "timing = " << timing; + + int totalpixels = 3712 * 3712; + + if(imageptrs->rgbrecipes[recipe].reflectivechannel) + { + time = new double[totalpixels]; + lat = new float[totalpixels]; + lon = new float[totalpixels]; + sza = new float[totalpixels]; + + snu_init_array_d(time, totalpixels, FILL_VALUE_F); + snu_init_array_f(lat, totalpixels, FILL_VALUE_F); + snu_init_array_f(lon, totalpixels, FILL_VALUE_F); + snu_init_array_f(sza, totalpixels, FILL_VALUE_F); + + saa = new float[totalpixels]; + vza = new float[totalpixels]; + vaa = new float[totalpixels]; + + snu_init_array_f(saa, totalpixels, FILL_VALUE_F); + snu_init_array_f(vza, totalpixels, FILL_VALUE_F); + snu_init_array_f(vaa, totalpixels, FILL_VALUE_F); + } + + result[0] = new float[totalpixels]; + result[1] = new float[totalpixels]; + result[2] = new float[totalpixels]; + + snu_init_array_f(result[0], totalpixels, FILL_VALUE_F); + snu_init_array_f(result[1], totalpixels, FILL_VALUE_F); + snu_init_array_f(result[2], totalpixels, FILL_VALUE_F); + + MsgFileAccess fa(directory, "H", productid1, productid2, timing); + MsgDataAccess da; + MSG_data prodata; + MSG_data epidata; + MSG_header header; + + QApplication::setOverrideCursor( Qt::WaitCursor ); // this might take time + + da.scan(fa, prodata, epidata, header); + emit progressCounter(20); + + double jtime_start = epidata.epilogue->product_stats.ActualScanningSummary.ForwardScanStart.get_jtime(); + double jtime_end = epidata.epilogue->product_stats.ActualScanningSummary.ForwardScanEnd.get_jtime(); + + double jtime = (jtime_start + jtime_end) / 2.; + + struct tm cdate; + Calendar_Date(jtime, &cdate); + double day_of_year = jtime - QSgp4Date::DateToJD(cdate.tm_year, 1, 0, true); + + qDebug() << "Julian day = " << jtime; + qDebug() << "day_of_year = " << day_of_year; + + double jtime2; + double jtime_start2; + double jtime_end2; + + double slope, offset; + bool toint; + + for(int i = 0; i < uniquebandlist.length(); i++) + { + + bandstorage newband; + newband.listindex = i; + newband.spectral_channel_nbr = uniquechannelnbrlist.at(i); + prodata.prologue->radiometric_proc.get_slope_offset(newband.spectral_channel_nbr, slope, offset, toint); + + newband.data = new float[totalpixels]; + snu_init_array_f(newband.data, totalpixels, FILL_VALUE_F); + newband.min = 0; newband.max = 0; + newband.directory = directory; + newband.productid1 = productid1; + newband.productid2 = uniquebandlist.at(i); + newband.timing = timing; + newband.day_of_year = day_of_year; + newband.units = uniqueunitlist.at(i); + newband.slope = slope; + newband.offset = offset; + bands.append(newband); + } + + if(imageptrs->rgbrecipes[recipe].reflectivechannel) + { + + /*------------------------------------------------------------------------- + * Compute the satellite position vector in Cartesian coordinates (km). + *-----------------------------------------------------------------------*/ + int i; + for (i = 0; i < 100; ++i) { + jtime_start2 = prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].StartTime.get_jtime(); + jtime_end2 = prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].EndTime.get_jtime(); + if (jtime >= jtime_start2 && jtime <= jtime_end2) + break; + } + + if (i == 100) { + fprintf(stderr, "ERROR: Image time is out of range of supplied orbit " + "polynomials\n"); + return; + } + + double t, X, Y, Z; + t = (jtime - (jtime_start2 + jtime_end2) / 2.) / ((jtime_end2 - jtime_start2) / 2.); + + X = prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].X[0] + + prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].X[1] * t; + Y = prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].Y[0] + + prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].Y[1] * t; + Z = prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].Z[0] + + prodata.prologue->sat_status.Orbit.OrbitPoliniomal[i].Z[1] * t; + /*------------------------------------------------------------------------- + * Compute latitude and longitude and solar and sensor zenith and azimuth + * angles. + *-----------------------------------------------------------------------*/ + + double lon0 = prodata.prologue->image_description.ProjectionDescription.LongitudeOfSSP; + qDebug() << "longitude SSP = " << lon0; + + double mu0; + double theta0; + double phi0; + + quint16 seglines = header.image_structure->number_of_lines; + quint16 columns = header.image_structure->number_of_columns; + + unsigned int i_image; + + for (int i = 0; i < 3712; ++i) + { + jtime2 = jtime_start + (double) i / (double) (3712 - 1) * (jtime_end - jtime_start); + for (int j = 0; j < 3712; ++j) + { + i_image = i * 3712 + j; + + snu_line_column_to_lat_lon(i, j, &lat[i_image], &lon[i_image], lon0, &nav_scaling_factors_vir); + + if (lat[i_image] != FILL_VALUE_F && lon[i_image] != FILL_VALUE_F) + { + time[i_image] = jtime2; + + snu_solar_params2(jtime2, lat[i_image] * D2R, lon[i_image] * D2R, &mu0, &theta0, &phi0, NULL); + sza[i_image] = theta0 * R2D; + saa[i_image] = phi0 * R2D; + + saa[i_image] = saa[i_image] + 180.; + if (saa[i_image] > 360.) + saa[i_image] = saa[i_image] - 360.; + + snu_vza_and_vaa(lat[i_image], lon[i_image], 0., X, Y, Z, &vza[i_image], &vaa[i_image]); + + vaa[i_image] = vaa[i_image] + 180.; + if (vaa[i_image] > 360.) + vaa[i_image] = vaa[i_image] - 360.; + + } + else + { + for(int k = 0; k < bands.length(); k++) + bands[k].data[i_image] = FILL_VALUE_F; + } + } + } + + PrintResults(sza, "sza"); + PrintResults(vza, "vza"); + } + + emit progressCounter(30); + + if(recipe != 4 && recipe != 6) + { + QtConcurrent::blockingMap(bands, [this] (bandstorage &data) { CalculateGeoRadiances(data); }); + + //Printbands(); + + bool pixelok = false; + + for(int colorindex = 0; colorindex < 3; colorindex++) + { + emit progressCounter(40 + colorindex * 10); + + for(int i = 0; i < imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).channels.length(); i++) + { + for(int j = 0; j < bands.length(); j++) + { + if(bands[j].spectral_channel_nbr == imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).spectral_channel_nbr.at(i)) + { + qDebug() << colorindex << " " << imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).channels.at(i) << " " << + imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).spectral_channel_nbr.at(i) << " " << + imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).subtract.at(i); + for (int line = 0; line < 3712; line++) + { + for (int pixelx = 0; pixelx < 3712; pixelx++) + { + int i_image = line * 3712 + pixelx; + pixelok = true; + for(int j1 = 0; j1 < bands.length();j1++) + { + if(bands[j1].data[i_image] == FILL_VALUE_F || bands[j1].data[i_image] == FILL_VALUE_US ) + { + pixelok = false; + break; + } + } + if(pixelok) + { + if(imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).subtract.at(i)) + { + if(result[colorindex][i_image] == FILL_VALUE_F) + result[colorindex][i_image] = - bands[j].data[i_image]; + else + result[colorindex][i_image] -= bands[j].data[i_image]; + } + else + { + if(result[colorindex][i_image] == FILL_VALUE_F) + result[colorindex][i_image] = bands[j].data[i_image]; + else + result[colorindex][i_image] += bands[j].data[i_image]; + } + } + else + { + result[colorindex][i_image] = FILL_VALUE_F; + } + } + } + } + } + } + } + } + else + { + if(recipe == 4) + { + ComposeDayMicrophysicsRGB(bands[0], jtime); + } + else if(recipe == 6) + { + ComposeSnowRGB(bands[0], jtime); + } + } + //PrintResults(); + + float resultmax[3]; + float resultmin[3]; + long countresultfill[3]; + + for(int i = 0; i < 3; i++) + { + resultmin[i] = 999999999; + resultmax[i] = -999999999; + countresultfill[i] = 0; + } + + emit progressCounter(70); + + for(int colorindex = 0; colorindex < 3; colorindex++) + { + for (int line = 0; line < 3712; line++) + { + for (int pixelx = 0; pixelx < 3712; pixelx++) + { + int i_image = line * 3712 + pixelx; + + if(result[colorindex][i_image] != FILL_VALUE_F) + { + if(result[colorindex][i_image] < resultmin[colorindex]) resultmin[colorindex] = result[colorindex][i_image]; + if(result[colorindex][i_image] > resultmax[colorindex]) resultmax[colorindex] = result[colorindex][i_image]; + } + else + countresultfill[colorindex]++; + } + } + } + + for(int colorindex = 0; colorindex < 3; colorindex++ ) + { + qDebug() << QString("---> %1 resultmin = %2 resultmax = %3").arg(colorindex).arg(resultmin[colorindex]).arg(resultmax[colorindex]); + qDebug() << QString("---> countresultfill = %1").arg(countresultfill[colorindex]); + } + + QRgb *row_col; + float red, green, blue; + imageptrs->InitializeImageGeostationary(opts.geosatellites.at(0).imagewidth, opts.geosatellites.at(0).imageheight); + + + for(int colorindex = 0; colorindex < 3; colorindex++) + { + for (int line = 0; line < 3712; line++) + { + for (int pixelx = 0; pixelx < 3712; pixelx++) + { + int i_image = line * 3712 + pixelx; + if(result[colorindex][i_image] != FILL_VALUE_F ) + { + if(imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).dimension == "K") + { + if(result[colorindex][i_image] > imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).rangeto ) + result[colorindex][i_image] = imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).rangeto; + if(result[colorindex][i_image] < imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).rangefrom ) + result[colorindex][i_image] = imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).rangefrom; + } + if(imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).dimension == "%") + { + if(result[colorindex][i_image] > resultmax[colorindex]*imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).rangeto) + result[colorindex][i_image] = resultmax[colorindex]*imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).rangeto; + + } + } + } + } + } + + for(int colorindex = 0; colorindex < 3; colorindex++ ) + { + if(imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).dimension == "K") + { + resultmin[colorindex] = imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).rangefrom; + resultmax[colorindex] = imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).rangeto; + } + if(imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).dimension == "%") + { + resultmax[colorindex] = resultmax[colorindex]*imageptrs->rgbrecipes[recipe].Colorvector.at(colorindex).rangeto; + } + } + + emit progressCounter(80); + + for(int colorindex = 0; colorindex < 3; colorindex++ ) + qDebug() << QString("%1 resultmin = %2 resultmax = %3").arg(colorindex).arg(resultmin[colorindex]).arg(resultmax[colorindex]); + + for (int line = 0; line < 3712; line++) + { + row_col = (QRgb*)imageptrs->ptrimageGeostationary->scanLine(3711 - line); + for (int pixelx = 0; pixelx < 3712; pixelx++) + { + int i_image = line * 3712 + pixelx; + + if(result[0][i_image] == FILL_VALUE_F || result[1][i_image] == FILL_VALUE_F || result[2][i_image] == FILL_VALUE_F ) + { + red = 0.0; + green = 0.0; + blue = 0.0; + } + else + { + if(imageptrs->rgbrecipes[recipe].Colorvector.at(0).inverse.at(0)) + red = 255.0 - 255.0 * pow((result[0][i_image] - resultmin[0]) / (resultmax[0] - resultmin[0]), 1.0/imageptrs->rgbrecipes[recipe].Colorvector.at(0).gamma); + else + red = 255.0 * pow((result[0][i_image] - resultmin[0]) / (resultmax[0] - resultmin[0]), 1.0/imageptrs->rgbrecipes[recipe].Colorvector.at(0).gamma); + if(imageptrs->rgbrecipes[recipe].Colorvector.at(1).inverse.at(0)) + green = 255.0 - 255.0 * pow((result[1][i_image] - resultmin[1]) / (resultmax[1] - resultmin[1]), 1.0/imageptrs->rgbrecipes[recipe].Colorvector.at(1).gamma); + else + green = 255.0 * pow((result[1][i_image] - resultmin[1]) / (resultmax[1] - resultmin[1]), 1.0/imageptrs->rgbrecipes[recipe].Colorvector.at(1).gamma); + if(imageptrs->rgbrecipes[recipe].Colorvector.at(2).inverse.at(0)) + blue = 255 - 255.0 * pow((result[2][i_image] - resultmin[2]) / (resultmax[2] - resultmin[2]), 1.0/imageptrs->rgbrecipes[recipe].Colorvector.at(2).gamma); + else + blue = 255.0 * pow((result[2][i_image] - resultmin[2]) / (resultmax[2] - resultmin[2]), 1.0/imageptrs->rgbrecipes[recipe].Colorvector.at(2).gamma); + } + + row_col[3711 - pixelx] = qRgb((int)red, (int)green, (int)blue); + } + } + + + qDebug() << "Kind of image = " << kindofimage; + emit signalcomposefinished(kindofimage); + emit progressCounter(100); + + ////////////////////////////////////////////////// + // Cleaup + ////////////////////////////////////////////////// + + for(int i = 0; i < uniquebandlist.length(); i++) + { + delete [] bands[i].data; + } + + bands.clear(); + + if(imageptrs->rgbrecipes[recipe].reflectivechannel) + { + delete [] time; + delete [] lat; + delete [] lon; + delete [] sza; + delete [] saa; + delete [] vza; + delete [] vaa; + } + + delete [] result[0]; + delete [] result[1]; + delete [] result[2]; + + qDebug() << "CreateRGBrecipeImageThread Finished !"; +} + +void SegmentListGeostationary::ComposeDayMicrophysicsRGB(bandstorage &bs, double julian_day) +{ + +// Red VIS008 0 --> 100 % Gamma 1.0 +// Green IR3.9Refl 0 --> 60 % Gamma 2.5 +// Blue IR10.8 203 --> 323 K Gamma 1.0 + + //VIS006 1 + //VIS008 2 + //IR_016 3 + //IR_039 4 + //WV_062 5 + //WV_073 6 + //IR_087 7 + //IR_097 8 + //IR_108 9 + //IR_120 10 + //IR_134 11 + + // REFL(IR3.9) = 100 * (Rtot - Rtherm)/(TOARAD - Rtherm) + // + int totalpixels = 3712 * 3712; + + float *ref008 = new float[totalpixels]; + float *rad039 = new float[totalpixels]; + float *rad039_corr = new float[totalpixels]; + float *bt108 = new float[totalpixels]; + float *bt134 = new float[totalpixels]; + float *r_therm = new float[totalpixels]; + float *toarad = new float[totalpixels]; + + snu_init_array_f(ref008, totalpixels, FILL_VALUE_F); + snu_init_array_f(rad039, totalpixels, FILL_VALUE_F); + snu_init_array_f(rad039_corr, totalpixels, FILL_VALUE_F); + snu_init_array_f(bt108, totalpixels, FILL_VALUE_F); + snu_init_array_f(bt134, totalpixels, FILL_VALUE_F); + snu_init_array_f(r_therm, totalpixels, FILL_VALUE_F); + snu_init_array_f(toarad, totalpixels, FILL_VALUE_F); + +//VIS008 + MsgFileAccess favis008(bs.directory, "H", bs.productid1, "VIS008", bs.timing); + MsgDataAccess da; + MSG_header header; + MSG_data prodata; + MSG_data epidata; + + qDebug() << "scanning vis008"; + + emit progressCounter(40); + +// da.scan(fac, header); + da.scan(favis008, prodata, epidata, header); + + double slope, offset; + bool toint; + + int channel = 2; // VIS008 + prodata.prologue->radiometric_proc.get_slope_offset(channel, slope, offset, toint); + for (size_t j = 0; j < da.segnames.size(); ++j) + { + cout << "Segment " << j << ": "; + MSG_data* d = da.segment(j); + if (!d) + { + cout << "missing.\n"; + continue; + } else { + for (size_t k = 0; k < da.npixperseg; ++k) + { + if(d->image->data[k] > 0) + { + ref008[j*da.npixperseg + k] = d->image->data[k]; + } + } + cout << std::flush; + } + } + + double a = PI / snu_solar_distance_factor2(bs.day_of_year); + int satid = (int)header.segment_id->spacecraft_id - 321 ; + double b = a / band_solar_irradiance[satid][channel - 1]; + + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (ref008[i_image] != FILL_VALUE_F && ref008[i_image] > 0) + { + if(sza[i_image] >= 0. && sza[i_image] < 80.) + { + double R = ref008[i_image] * slope + offset; + ref008[i_image] = b * R; + ref008[i_image] = ref008[i_image] < 0 ? 0 : ref008[i_image]; + } + else + ref008[i_image] = FILL_VALUE_F; + } + } + } + + emit progressCounter(45); + + //PrintResults(ref008, "ref008"); + + qDebug() << "scanning ir_108"; + + //BT from 10.8 + MsgFileAccess fair_108(bs.directory, "H", bs.productid1, "IR_108", bs.timing); + da.scan(fair_108, header); + + channel = 9; // IR_108 + prodata.prologue->radiometric_proc.get_slope_offset(channel, slope, offset, toint); + for (size_t j = 0; j < da.segnames.size(); ++j) + { + cout << "Segment " << j << ": "; + MSG_data* d = da.segment(j); + if (!d) + { + cout << "missing.\n"; + continue; + } else { + for (size_t k = 0; k < da.npixperseg; ++k) + { + if(d->image->data[k] > 0) + { + bt108[j*da.npixperseg + k] = d->image->data[k]; + } + } + cout << std::flush; + + } + } + + const double c1 = 1.19104e-5; // [mW m^-2 sr^-1 (cm^-1)^-4 + const double c2 = 1.43877; // K (cm^-1)^-1 + double nu = bt_nu_c[satid][channel - 1]; + double nu3 = nu * nu * nu; + + a = bt_A[satid][channel - 1]; + b = bt_B[satid][channel - 1]; + + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (bt108[i_image] != FILL_VALUE_F && bt108[i_image] > 0) + { + double L = bt108[i_image] * slope + offset; + + bt108[i_image] = (c2 * nu / log(1. + nu3 * c1 / L) - b) / a; + } + } + } + + emit progressCounter(50); + + // PrintResults(bt108, "bt108"); + + qDebug() << "scanning ir_134"; + + //BT from 13.4 + MsgFileAccess fair_134(bs.directory, "H", bs.productid1, "IR_134", bs.timing); + da.scan(fair_134, header); + + channel = 11; // IR_134 + prodata.prologue->radiometric_proc.get_slope_offset(channel, slope, offset, toint); + for (size_t j = 0; j < da.segnames.size(); ++j) + { + cout << "Segment " << j << ": "; + MSG_data* d = da.segment(j); + if (!d) + { + cout << "missing.\n"; + continue; + } else { + for (size_t k = 0; k < da.npixperseg; ++k) + { + if(d->image->data[k] > 0) + { + bt134[j*da.npixperseg + k] = d->image->data[k]; + } + } + cout << std::flush; + + } + } + + + nu = bt_nu_c[satid][channel - 1]; + nu3 = nu * nu * nu; + + a = bt_A[satid][channel - 1]; + b = bt_B[satid][channel - 1]; + + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (bt108[i_image] != FILL_VALUE_F && bt108[i_image] > 0) + { + double L = bt134[i_image] * slope + offset; + bt134[i_image] = (c2 * nu / log(1. + nu3 * c1 / L) - b) / a; + } + } + } + + //PrintResults(bt134, "bt134"); + + emit progressCounter(55); + + qDebug() << "calculating radiance ir_039"; + + // Radiance from IR_039 channel = 4 + // + channel = 4; + a = bt_A[satid][channel - 1]; + b = bt_B[satid][channel - 1]; + nu = bt_nu_c[satid][channel - 1]; + nu3 = nu * nu * nu; + + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (bt108[i_image] != FILL_VALUE_F && bt108[i_image] > 0) + { + rad039[i_image] = c1 * nu3 / (exp( c2 * nu / (a * bt108[i_image] + b)) - 1); + } + } + } + + emit progressCounter(60); + + //PrintResults(rad039, "calculating radiance ir_039"); + + + qDebug() << "calculating co2 correction"; + + //CO2 correction + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (bt108[i_image] != FILL_VALUE_F && bt108[i_image] > 0) + { + rad039_corr[i_image] = pow((bt108[i_image] - 0.25 * (bt108[i_image] - bt134[i_image])), 4.0) / pow(bt108[i_image], 4.0); + } + } + } + + //PrintResults(rad039_corr, "calculating co2 correction"); + + qDebug() << "calculating r_therm"; + + //r_therm = rad039 * rad039_corr + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (bt108[i_image] != FILL_VALUE_F && bt108[i_image] > 0) + { + r_therm[i_image] = rad039[i_image] * rad039_corr[i_image]; + } + } + } + + //PrintResults(r_therm, "calculating r_therm"); + + qDebug() << "calculating co2 correction solar constant"; + + //The CO2-corrected, solar constant at the Top of the Atmosphere + //in Channel 04 (IR3.9) + double ESD = 1.0 - 0.0167 * cos( 2 * PI * (julian_day - 3)/ 365.0); + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + if (rad039_corr[i_image] != FILL_VALUE_F && rad039_corr[i_image] > 0 && sza[i_image] < 80.0) + { + toarad[i_image] = 4.92 / ESD * ESD; + float costeta = cos(sza[i_image] * D2R); + float cossat = cos(vza[i_image] * D2R); + toarad[i_image] *= costeta; + + toarad[i_image] *= exp(-(1.0 - rad039_corr[i_image])); + toarad[i_image] *= exp(-(1.0 - rad039_corr[i_image]) * costeta / cossat); + } + } + } + + //PrintResults(toarad, "calculating co2 correction solar constant"); + + emit progressCounter(65); + + qDebug() << "calculating result"; + + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (rad039[i_image] != FILL_VALUE_F && rad039[i_image] > 0) + { + result[0][i_image] = ref008[i_image]; + //result[1][i_image] = 100 * (rad039[i_image] - r_therm[i_image]) / (toarad[i_image] - r_therm[i_image]); + result[1][i_image] = 100 * (rad039[i_image] - r_therm[i_image]); + result[2][i_image] = bt108[i_image]; + } + } + } + + //PrintResults(result[1], "green result"); + + delete [] ref008; + delete [] rad039; + delete [] rad039_corr; + delete [] bt108; + delete [] bt134; + delete [] r_therm; + delete [] toarad; + +} + +void SegmentListGeostationary::ComposeSnowRGB(bandstorage &bs, double julian_day) +{ + // Red VIS008 0 --> 100 % Gamma 1.7 + // Green IR_016 0 --> 70 % Gamma 1.7 + // Blue IR3.9Refl 0 --> 30 % Gamma 1.7 + + int totalpixels = 3712 * 3712; + + float *ref008 = new float[totalpixels]; + float *ref016 = new float[totalpixels]; + float *rad039 = new float[totalpixels]; + float *rad039_corr = new float[totalpixels]; + float *bt108 = new float[totalpixels]; + float *bt134 = new float[totalpixels]; + float *r_therm = new float[totalpixels]; + float *toarad = new float[totalpixels]; + + snu_init_array_f(ref008, totalpixels, FILL_VALUE_F); + snu_init_array_f(ref016, totalpixels, FILL_VALUE_F); + snu_init_array_f(rad039, totalpixels, FILL_VALUE_F); + snu_init_array_f(rad039_corr, totalpixels, FILL_VALUE_F); + snu_init_array_f(bt108, totalpixels, FILL_VALUE_F); + snu_init_array_f(bt134, totalpixels, FILL_VALUE_F); + snu_init_array_f(r_therm, totalpixels, FILL_VALUE_F); + snu_init_array_f(toarad, totalpixels, FILL_VALUE_F); + +//VIS008 + MsgFileAccess favis008(bs.directory, "H", bs.productid1, "VIS008", bs.timing); + MsgDataAccess da; + MSG_header header; + MSG_data prodata; + MSG_data epidata; + + qDebug() << "scanning vis008"; + + emit progressCounter(40); + +// da.scan(fac, header); + da.scan(favis008, prodata, epidata, header); + + double slope, offset; + bool toint; + + int channel = 2; // VIS008 + prodata.prologue->radiometric_proc.get_slope_offset(channel, slope, offset, toint); + for (size_t j = 0; j < da.segnames.size(); ++j) + { + cout << "Segment " << j << ": "; + MSG_data* d = da.segment(j); + if (!d) + { + cout << "missing.\n"; + continue; + } else { + for (size_t k = 0; k < da.npixperseg; ++k) + { + if(d->image->data[k] > 0) + { + ref008[j*da.npixperseg + k] = d->image->data[k]; + } + } + cout << std::flush; + } + } + + double a = PI / snu_solar_distance_factor2(bs.day_of_year); + int satid = (int)header.segment_id->spacecraft_id - 321 ; + double b = a / band_solar_irradiance[satid][channel - 1]; + + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (ref008[i_image] != FILL_VALUE_F && ref008[i_image] > 0) + { + if(sza[i_image] >= 0. && sza[i_image] < 80.) + { + double R = ref008[i_image] * slope + offset; + ref008[i_image] = b * R; + ref008[i_image] = ref008[i_image] < 0 ? 0 : ref008[i_image]; + } + else + ref008[i_image] = FILL_VALUE_F; + } + } + } + + //IR_016 + MsgFileAccess fair_016(bs.directory, "H", bs.productid1, "IR_016", bs.timing); + + qDebug() << "scanning ir_016"; + + da.scan(fair_016, header); + + channel = 3; // IR_016 + prodata.prologue->radiometric_proc.get_slope_offset(channel, slope, offset, toint); + for (size_t j = 0; j < da.segnames.size(); ++j) + { + cout << "Segment " << j << ": "; + MSG_data* d = da.segment(j); + if (!d) + { + cout << "missing.\n"; + continue; + } else { + for (size_t k = 0; k < da.npixperseg; ++k) + { + if(d->image->data[k] > 0) + { + ref016[j*da.npixperseg + k] = d->image->data[k]; + } + } + cout << std::flush; + } + } + + b = a / band_solar_irradiance[satid][channel - 1]; + + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (ref016[i_image] != FILL_VALUE_F && ref016[i_image] > 0) + { + if(sza[i_image] >= 0. && sza[i_image] < 80.) + { + double R = ref016[i_image] * slope + offset; + ref016[i_image] = b * R; + ref016[i_image] = ref016[i_image] < 0 ? 0 : ref016[i_image]; + } + else + ref016[i_image] = FILL_VALUE_F; + } + } + } + + emit progressCounter(45); + + //PrintResults(ref016, "ref016"); + + qDebug() << "scanning ir_108"; + + //BT from 10.8 + MsgFileAccess fair_108(bs.directory, "H", bs.productid1, "IR_108", bs.timing); + da.scan(fair_108, header); + + channel = 9; // IR_108 + prodata.prologue->radiometric_proc.get_slope_offset(channel, slope, offset, toint); + for (size_t j = 0; j < da.segnames.size(); ++j) + { + cout << "Segment " << j << ": "; + MSG_data* d = da.segment(j); + if (!d) + { + cout << "missing.\n"; + continue; + } else { + for (size_t k = 0; k < da.npixperseg; ++k) + { + if(d->image->data[k] > 0) + { + bt108[j*da.npixperseg + k] = d->image->data[k]; + } + } + cout << std::flush; + + } + } + + const double c1 = 1.19104e-5; // [mW m^-2 sr^-1 (cm^-1)^-4 + const double c2 = 1.43877; // K (cm^-1)^-1 + + double nu = bt_nu_c[satid][channel - 1]; + double nu3 = nu * nu * nu; + + a = bt_A[satid][channel - 1]; + b = bt_B[satid][channel - 1]; + + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (bt108[i_image] != FILL_VALUE_F && bt108[i_image] > 0) + { + double L = bt108[i_image] * slope + offset; + + bt108[i_image] = (c2 * nu / log(1. + nu3 * c1 / L) - b) / a; + } + } + } + + emit progressCounter(50); + + // PrintResults(bt108, "bt108"); + + qDebug() << "scanning ir_134"; + + //BT from 13.4 + MsgFileAccess fair_134(bs.directory, "H", bs.productid1, "IR_134", bs.timing); + da.scan(fair_134, header); + + channel = 11; // IR_134 + prodata.prologue->radiometric_proc.get_slope_offset(channel, slope, offset, toint); + for (size_t j = 0; j < da.segnames.size(); ++j) + { + cout << "Segment " << j << ": "; + MSG_data* d = da.segment(j); + if (!d) + { + cout << "missing.\n"; + continue; + } else { + for (size_t k = 0; k < da.npixperseg; ++k) + { + if(d->image->data[k] > 0) + { + bt134[j*da.npixperseg + k] = d->image->data[k]; + } + } + cout << std::flush; + + } + } + + + nu = bt_nu_c[satid][channel - 1]; + nu3 = nu * nu * nu; + + a = bt_A[satid][channel - 1]; + b = bt_B[satid][channel - 1]; + + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (bt108[i_image] != FILL_VALUE_F && bt108[i_image] > 0) + { + double L = bt134[i_image] * slope + offset; + + bt134[i_image] = (c2 * nu / log(1. + nu3 * c1 / L) - b) / a; + } + } + } + + //PrintResults(bt134, "bt134"); + + emit progressCounter(55); + + qDebug() << "calculating radiance ir_039"; + + // Radiance from IR_039 channel = 4 + // + channel = 4; + a = bt_A[satid][channel - 1]; + b = bt_B[satid][channel - 1]; + nu = bt_nu_c[satid][channel - 1]; + nu3 = nu * nu * nu; + + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (bt108[i_image] != FILL_VALUE_F && bt108[i_image] > 0) + { + rad039[i_image] = c1 * nu3 / (exp( c2 * nu / (a * bt108[i_image] + b)) - 1); + } + } + } + + emit progressCounter(60); + + //PrintResults(rad039, "calculating radiance ir_039"); + + + qDebug() << "calculating co2 correction"; + + //CO2 correction + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (bt108[i_image] != FILL_VALUE_F && bt108[i_image] > 0) + { + rad039_corr[i_image] = pow((bt108[i_image] - 0.25 * (bt108[i_image] - bt134[i_image])), 4.0) / pow(bt108[i_image], 4.0); + } + } + } + + //PrintResults(rad039_corr, "calculating co2 correction"); + + qDebug() << "calculating r_therm"; + + //r_therm = rad039 * rad039_corr + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (bt108[i_image] != FILL_VALUE_F && bt108[i_image] > 0) + { + r_therm[i_image] = rad039[i_image] * rad039_corr[i_image]; + } + } + } + + //PrintResults(r_therm, "calculating r_therm"); + + qDebug() << "calculating co2 correction solar constant"; + + //The CO2-corrected, solar constant at the Top of the Atmosphere + //in Channel 04 (IR3.9) + double ESD = 1.0 - 0.0167 * cos( 2 * PI * (julian_day - 3)/ 365.0); + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + if (rad039_corr[i_image] != FILL_VALUE_F && rad039_corr[i_image] > 0 && sza[i_image] < 80.0) + { + toarad[i_image] = 4.92 / ESD * ESD; + float costeta = cos(sza[i_image] * D2R); + float cossat = cos(vza[i_image] * D2R); + toarad[i_image] *= costeta; + + toarad[i_image] *= exp(-(1.0 - rad039_corr[i_image])); + toarad[i_image] *= exp(-(1.0 - rad039_corr[i_image]) * costeta / cossat); + } + } + } + + //PrintResults(toarad, "calculating co2 correction solar constant"); + + emit progressCounter(65); + + qDebug() << "calculating result"; + + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (rad039[i_image] != FILL_VALUE_F && rad039[i_image] > 0) + { + result[0][i_image] = ref008[i_image]; + result[1][i_image] = ref016[i_image]; + result[2][i_image] = 100 * (rad039[i_image] - r_therm[i_image]); + } + } + } + + //PrintResults(result[1], "green result"); + + delete [] ref008; + delete [] ref016; + delete [] rad039; + delete [] rad039_corr; + delete [] bt108; + delete [] bt134; + delete [] r_therm; + delete [] toarad; +} + +void SegmentListGeostationary::Printbands() +{ + for(int k = 0; k < bands.length(); k++) + { + cout << "band " << k << " " << bands[k].productid2.toStdString() << endl; + for (int i = 1000; i < 1000 + 10; ++i) + { + for (int j = 1856; j < 1856 + 10; ++j) + { + int i_image = i * 3712 + j; + + cout << bands[k].data[i_image] << " "; + } + cout << endl; + } + cout << endl; + } + +// cout << "sza " << endl; +// for (int i = 1000; i < 1000 + 10; ++i) +// { +// for (int j = 1856; j < 1856 + 10; ++j) +// { +// int i_image = i * 3712 + j; + +// cout << sza[i_image] << " "; +// } +// cout << endl; +// } + +} + +void SegmentListGeostationary::PrintResults() +{ + for(int k = 0; k < 3; k++) + { + cout << "color " << k << endl; + for (int i = 1000; i < 1000 + 10; ++i) + { + for (int j = 1856; j < 1856 + 10; ++j) + { + int i_image = i * 3712 + j; + + cout << result[k][i_image] << " "; + } + cout << endl; + } + cout << endl; + } +} + +void SegmentListGeostationary::PrintResults(float *ptr, QString title) +{ + + cout << title.toStdString() << endl; + for (int i = 1000; i < 1000 + 10; ++i) + { + for (int j = 1856; j < 1856 + 10; ++j) + { + int i_image = i * 3712 + j; + + cout << ptr[i_image] << " "; + } + cout << endl; + } + cout << endl << std::flush; + +} + +void SegmentListGeostationary::CalculateGeoRadiances(bandstorage &bs) +{ + qDebug() << QString("ComposeGeoRGBRecipe %1 %2").arg(bs.listindex).arg(bs.productid2); + + MsgFileAccess fac(bs.directory, "H", bs.productid1, bs.productid2, bs.timing); + MsgDataAccess da; + MSG_header header; + + da.scan(fac, header); + MSG_SAMPLE compmin = 0xffff, compmax = 0; + for (size_t j = 0; j < da.segnames.size(); ++j) + { + cout << "Segment " << j << ": "; + MSG_data* d = da.segment(j); + MSG_SAMPLE min = 0xffff, max = 0; + if (!d) + { + cout << "missing.\n"; + continue; + } else { + for (size_t k = 0; k < da.npixperseg; ++k) + { + if(d->image->data[k] > 0) + { + bands[bs.listindex].data[j*da.npixperseg + k] = d->image->data[k]; + if (d->image->data[k] < min) min = d->image->data[k]; + if (d->image->data[k] > max) max = d->image->data[k]; + } + } + cout << "min " << min << " max " << max << endl; + } + if (min < compmin) compmin = min; + if (max > compmax) compmax = max; + } + cout << "compmin = " << compmin << " compmax = " << compmax << endl; + bands[bs.listindex].min = (float)compmin; + bands[bs.listindex].max = (float)compmax; + + + /*------------------------------------------------------------------------- + * Compute radiance for the bands requested. + * + * Ref: PDF_TEN_05105_MSG_IMG_DATA, Page 26 + *-----------------------------------------------------------------------*/ + + double slope, offset; + bool toint; + + if(bands[bs.listindex].units == SEVIRI_UNIT_RAD) // mW*pow(m,-2)*pow(sr,-1)*pow(pow(cm,-1)), -1) + { + bands[bs.listindex].min = std::numeric_limits::max(); + bands[bs.listindex].max = std::numeric_limits::min(); + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (bands[bs.listindex].data[i_image] != FILL_VALUE_US && bands[bs.listindex].data[i_image] > 0) + { + bands[bs.listindex].data[i_image] = bands[bs.listindex].data[i_image] * bs.slope + bs.offset; + bands[bs.listindex].data[i_image] = bands[bs.listindex].data[i_image] < 0 ? 0 : bands[bs.listindex].data[i_image]; + if(bands[bs.listindex].data[i_image] < bands[bs.listindex].min) bands[bs.listindex].min = bands[bs.listindex].data[i_image]; + if(bands[bs.listindex].data[i_image] > bands[bs.listindex].max) bands[bs.listindex].max = bands[bs.listindex].data[i_image]; + } + } + } + } + + + if(bands[bs.listindex].units == SEVIRI_UNIT_RAD) + { + qDebug() << QString("Radiances %1 minimum = %2 maximum = %3").arg(bands[bs.listindex].productid2).arg(bands[bs.listindex].min).arg(bands[bs.listindex].max); + } + + /*------------------------------------------------------------------------- + * Compute reflectance or bidirectional reflectance factor (BRF) for the + * bands requested. + * + * Ref: PDF_MSG_SEVIRI_RAD2REFL, Page 8 + *-----------------------------------------------------------------------*/ + //double dd = 1. / sqrt(snu_solar_distance_factor2(day_of_year)); // distance earth-sun + //double a = PI * dd * dd; + double a = PI / snu_solar_distance_factor2(bs.day_of_year); + + int satid = (int)header.segment_id->spacecraft_id - 321 ; + + if (bands[bs.listindex].units == SEVIRI_UNIT_REF || bands[bs.listindex].units == SEVIRI_UNIT_BRF) + { + bands[bs.listindex].min = std::numeric_limits::max(); + bands[bs.listindex].max = std::numeric_limits::min(); + + double b = a / band_solar_irradiance[satid][bands[bs.listindex].spectral_channel_nbr - 1]; + + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (bands[bs.listindex].units == SEVIRI_UNIT_BRF && bands[bs.listindex].data[i_image] != FILL_VALUE_US && bands[bs.listindex].data[i_image] > 0) + { + if(sza[i_image] >= 0. && sza[i_image] < 80.) + { + double R = bands[bs.listindex].data[i_image] * bs.slope + bs.offset; + + bands[bs.listindex].data[i_image] = b * R; + + if (bands[bs.listindex].units == SEVIRI_UNIT_BRF) + bands[bs.listindex].data[i_image] /= cos(sza[i_image] * D2R); + + bands[bs.listindex].data[i_image] = bands[bs.listindex].data[i_image] < 0.0 ? 0.0 : bands[bs.listindex].data[i_image]; + if( i_image == 1856 * 3712 + 1856 ) + qDebug() << QString("The effective radiance = %1 BRF = %2").arg(R).arg(bands[bs.listindex].data[i_image]); + + + if(bands[bs.listindex].data[i_image] < bands[bs.listindex].min) bands[bs.listindex].min = bands[bs.listindex].data[i_image]; + if(bands[bs.listindex].data[i_image] > bands[bs.listindex].max) bands[bs.listindex].max = bands[bs.listindex].data[i_image]; + } + else + bands[bs.listindex].data[i_image] = FILL_VALUE_US; + } + else if (bands[bs.listindex].units == SEVIRI_UNIT_REF && bands[bs.listindex].data[i_image] != FILL_VALUE_US && bands[bs.listindex].data[i_image] > 0) + { + double R = bands[bs.listindex].data[i_image] * bs.slope + bs.offset; + + bands[bs.listindex].data[i_image] = b * R; + + if(bands[bs.listindex].data[i_image] < bands[bs.listindex].min) bands[bs.listindex].min = bands[bs.listindex].data[i_image]; + if(bands[bs.listindex].data[i_image] > bands[bs.listindex].max) bands[bs.listindex].max = bands[bs.listindex].data[i_image]; + } + } + } + } + + const double c1 = 1.19104e-5; // [mW m^-2 sr^-1 (cm^-1)^-4 + const double c2 = 1.43877; // K (cm^-1)^-1 + + /*------------------------------------------------------------------------- + * Compute brightness temperature for the bands requested. + * + * Ref: PDF_TEN_05105_MSG_IMG_DATA, Page 26 + * Ref: The Conversion from Effective Radiances + * to Equivalent Brightness Temperatures (EUM/MET/TEN/11/0569) + *-----------------------------------------------------------------------*/ + if (bands[bs.listindex].units == SEVIRI_UNIT_BT) + { + bands[bs.listindex].min = std::numeric_limits::max(); + bands[bs.listindex].max = std::numeric_limits::min(); + + qDebug() << bands[bs.listindex].spectral_channel_nbr << " " << bands[bs.listindex].productid2 << " slope = " << bs.slope << " offset = " << bs.offset; + /* + nu = 1.e4 / channel_center_wavelength[d->image.band_ids[bs.listindex] - 1]; +*/ + double nu = bt_nu_c[satid][bands[bs.listindex].spectral_channel_nbr - 1]; + double nu3 = nu * nu * nu; + + double a = bt_A[satid][bands[bs.listindex].spectral_channel_nbr - 1]; + double b = bt_B[satid][bands[bs.listindex].spectral_channel_nbr - 1]; + + for (int j = 0; j < 3712; ++j) + { + for (int k = 0; k < 3712; ++k) + { + int i_image = j * 3712 + k; + + if (bands[bs.listindex].data[i_image] != FILL_VALUE_F && bands[bs.listindex].data[i_image] > 0) + { + double L = bands[bs.listindex].data[i_image] * bs.slope + bs.offset; + + bands[bs.listindex].data[i_image] = (c2 * nu / log(1. + nu3 * c1 / L) - b) / a; + if( i_image == 1856 * 3712 + 1856 ) + qDebug() << QString("The effective radiance = %1 the brightense temp = %2").arg(L).arg(bands[bs.listindex].data[i_image]); + + if(bands[bs.listindex].data[i_image] < bands[bs.listindex].min) bands[bs.listindex].min = bands[bs.listindex].data[i_image]; + if(bands[bs.listindex].data[i_image] > bands[bs.listindex].max) bands[bs.listindex].max = bands[bs.listindex].data[i_image]; + } + } + } + } + + qDebug() << QString("%1 minimum = %2 maximum = %3").arg(bands[bs.listindex].productid2).arg(bands[bs.listindex].min).arg(bands[bs.listindex].max); + +} diff --git a/core/segmentlistgeostationary.h b/core/segmentlistgeostationary.h index 40bf17b..e1a0efa 100644 --- a/core/segmentlistgeostationary.h +++ b/core/segmentlistgeostationary.h @@ -6,6 +6,26 @@ #include #include #include "globals.h" +#include "msgfileaccess.h" +#include "msgdataaccess.h" +#include "segmentimage.h" + +typedef struct { + int listindex; + int spectral_channel_nbr; + QString directory; + QString productid1; + QString productid2; + QString timing; + double day_of_year; + float min; + float max; + float *data; + int units; + double slope; + double offset; +} bandstorage; + class SegmentListGeostationary : public QObject { @@ -37,12 +57,12 @@ class SegmentListGeostationary : public QObject void ComposeSegmentImageHDF(QFileInfo fileinfo, int channelindex, QVector spectrumvector, QVector inversevector ); void ComposeSegmentImageHDFInThread(QStringList filelist, QVector spectrumvector, QVector inversevector ); - void ComposeSegmentImagenetCDFInThread(); //QStringList filelist, QVector spectrumvector, QVector inversevector, int histogrammethod ); - void SetupContrastStretch(quint16 x1, quint16 y1, quint16 x2, quint16 y2); //, quint16 x3, quint16 y3, quint16 x4, quint16 y4); + void ComposeSegmentImagenetCDFInThread(); + void SetupContrastStretch(quint16 x1, quint16 y1, quint16 x2, quint16 y2); quint16 ContrastStretch(quint16 val); void InsertPresent( QVector spectrumvector, QString filespectrum, int filesequence); bool allHRVColorSegmentsReceived(); - bool allSegmentsReceived(); + bool allVIS_IRSegmentsReceived(); bool bActiveSegmentList; bool bisRSS; eGeoSatellite getGeoSatellite(); @@ -52,11 +72,17 @@ class SegmentListGeostationary : public QObject void CalculateLUTGeo(int colorindex); void CalculateLUTGeo(int colorindex, quint16 *ptr, quint16 fillvalue); + void ComposeGeoRGBRecipe(int recipe, QString tex); + void ComposeGeoRGBRecipeInThread(int recipe); static void doComposeGeostationaryXRIT(SegmentListGeostationary *sm, QString segment_path, int channelindex, QVector spectrumvector, QVector inversevector); static void doComposeGeostationaryXRITHimawari(SegmentListGeostationary *sm, QString segment_path, int channelindex, QVector spectrumvector, QVector inversevector); static void doComposeGeostationaryHDFInThread(SegmentListGeostationary *sm, QStringList filelist, QVector spectrumvector, QVector inversevector); - static void doComposeGeostationarynetCDFInThread(SegmentListGeostationary *sm); //, QStringList strlist, QVector spectrumvector, QVector inversevector, int histogrammethod); + static void doComposeGeostationarynetCDFInThread(SegmentListGeostationary *sm); + static void doComposeGeoRGBRecipe(SegmentListGeostationary *sm, int recipe); + void CalculateGeoRadiances(bandstorage &bs); void setThreadParameters(QStringList strlist, QVector spectrumvector, QVector inversevector, int histogrammethod, bool pseudocolor); + void ComposeDayMicrophysicsRGB(bandstorage &bs, double julian_day); + void ComposeSnowRGB(bandstorage &bs, double julian_day); bool issegmentcomposedRed[10]; bool issegmentcomposedGreen[10]; @@ -97,6 +123,9 @@ class SegmentListGeostationary : public QObject void ComposeVISIR(); void ComposeVISIRHimawari(); void getFilenameParameters(QFileInfo fileinfo, QString &filespectrum, QString &filedate, int &filesequence); + void Printbands(); + void PrintResults(); + void PrintResults(float *ptr, QString title); quint16 stat_min[3]; quint16 stat_max[3]; @@ -120,6 +149,17 @@ class SegmentListGeostationary : public QObject int histogrammethod; bool pseudocolor; + QList bands; + float* result[3]; + double *time; + float *lat; + float *lon; + + float *sza; /* image of solar zenith angle (degrees: 0.0 -- 180.0) */ + float *saa; /* image of solar azimuth angle (degrees: 0.0 -- 360.0) */ + float *vza; /* image of viewing zenith angle (degrees: 0.0 -- 180.0) */ + float *vaa; /* image of viewing azimuth angle (degrees: 0.0 -- 360.0) */ + QString tex; diff --git a/core/stereographic.cpp b/core/stereographic.cpp index 47e008a..2ba2095 100644 --- a/core/stereographic.cpp +++ b/core/stereographic.cpp @@ -205,7 +205,7 @@ void StereoGraphic::CreateMapFromGeostationary() } else { - if(sl->getGeoSatellite() == eGeoSatellite::MET_9 || sl->getGeoSatellite() == eGeoSatellite::MET_10) + if(sl->getGeoSatellite() == eGeoSatellite::MET_10) { if( picrow >= 0 && picrow < 5*464) {