OSGeo
MSGN: avoid 'Scanline corrupt' report for HRV with or without RSS
Using archive data from https://data.eumetsat.int/data/map/EO:EUM:DAT:MSG:HRSEVIRI, whole earth images are prefixed MSG4-, RSS (every 5 minutes, only Europe) are MSG3-
To extract the HRV image, filename must be given as HRV:
Commit 84e31284211509314c84a41d78bd12962b9b170d allowed extraction of HRV images from whole-disk (MSG4-) files
RSS HRV (MSG3-) reports 'Scanline corrupt' - in fact the extracted image would be OK
This commit works around poGDS->m_bHRVDealWithSplit being set false for HRV RSS images (in msgndataset.cpp).
What does this PR do?
What are related issues/pull requests?
Tasklist
- [ ] ADD YOUR TASKS HERE
- [ ] Add test case(s)
- [ ] Add documentation
- [ ] Review
- [ ] Adjust for comments
- [ ] All CI builds and checks have passed
Environment
Provide environment details, if relevant:
- OS:
- Compiler:
Dear Even
Thank you for improving what I wrote - I wasn't aware of the naming convention.
As you can tell, this is my first time contributing to GDAL - git is familiar to me, but all of the github machinery is new to me. Tell me if there are other things I need to do to to finish this pull request.
Thanks again
--David
On Sun, 05 Jun 2022 02:07:35 -0700 Even Rouault @.***> wrote:
@rouault commented on this pull request.
@@ -61,16 +67,16 @@ class MSGNDataset final: public GDALDataset
Msg_reader_core* msg_reader_core; open_mode_type m_open_mode = MODE_VISIR;
- bool m_bHRVDealWithSplit = false;
- hrv_mode_type m_bHRVDealWithSplit = NOT_HRV;
The 'b' means boolean, hence changing to 'e' to mean enumeration:
hrv_mode_type m_eHRVDealWithSplit = NOT_HRV;@@ -209,7 +215,7 @@ CPLErr MSGNRasterBand::IReadBlock( CPL_UNUSED int nBlockXOff, SUB_VISIRLINE* p = (SUB_VISIRLINE*) pszRecord; to_native(*p);
- if (p->lineValidity != 1 || poGDS->m_bHRVDealWithSplit) {
- if (p->lineValidity != 1 || poGDS->m_bHRVDealWithSplit) { // Split lines are not full width, so NODATA all first
if (p->lineValidity != 1 || poGDS->m_eHRVDealWithSplit != NOT_HRV) { // Split lines are not full width, so NODATA all first@@ -220,8 +226,19 @@ CPLErr MSGNRasterBand::IReadBlock( CPL_UNUSED int nBlockXOff, }
if ( nread != data_length ||
(p->lineNumberInVisirGrid - poGDS->msg_reader_core->get_line_start()) != (unsigned int)i_nBlockYOff )
(p->lineNumberInVisirGrid - ((open_mode == MODE_HRV && poGDS->m_bHRVDealWithSplit == HRV_RSS) ?(p->lineNumberInVisirGrid - ((open_mode == MODE_HRV && poGDS->m_eHRVDealWithSplit == HRV_RSS) ?{
CPLDebug("MSGN", "Deal with split %s", poGDS->m_bHRVDealWithSplit == HRV_RSS ? "RSS" :CPLDebug("MSGN", "Deal with split %s", poGDS->m_eHRVDealWithSplit == HRV_RSS ? "RSS" :{
CPLDebug("MSGN", "Deal with split %s", poGDS->m_bHRVDealWithSplit == HRV_RSS ? "RSS" :(poGDS->m_bHRVDealWithSplit ? "whole" : "no") );(poGDS->m_eHRVDealWithSplit ? "whole" : "no") );@@ -472,17 +494,40 @@ GDALDataset *MSGNDataset::Open( GDALOpenInfo * poOpenInfo ) idr.plannedCoverage_hrv.upperWestColumnPlanned <= poDS->nRasterXSize * 3 ) { poDS->nRasterXSize *= 3;
poDS->m_bHRVDealWithSplit = true;
poDS->m_bHRVDealWithSplit = HRV_WHOLE_DISK;poDS->m_eHRVDealWithSplit = HRV_WHOLE_DISK;
idr.plannedCoverage_hrv.upperSouthLinePlanned == 0 &&
idr.plannedCoverage_hrv.upperWestColumnPlanned == 0 &&idr.plannedCoverage_hrv.upperEastColumnPlanned == 0 && // RSS only uses the lower sectionidr.plannedCoverage_hrv.lowerNorthLinePlanned == idr.referencegrid_hrv.numberOfLines && // start at max N// full expected widthidr.plannedCoverage_hrv.lowerWestColumnPlanned == idr.plannedCoverage_hrv.lowerEastColumnPlanned + nRawHRVColumns - 1 &&idr.plannedCoverage_hrv.lowerSouthLinePlanned > 1 &&idr.plannedCoverage_hrv.lowerSouthLinePlanned < idr.referencegrid_hrv.numberOfLines &&idr.plannedCoverage_hrv.lowerEastColumnPlanned >= 1 &&idr.plannedCoverage_hrv.lowerWestColumnPlanned <= poDS->nRasterXSize * 3 &&// full heightidr.plannedCoverage_hrv.lowerNorthLinePlanned == idr.plannedCoverage_hrv.lowerSouthLinePlanned + poDS->nRasterYSize -1){poDS->nRasterXSize *= 3;poDS->m_bHRVDealWithSplit = HRV_RSS;poDS->m_eHRVDealWithSplit = HRV_RSS;else{ // In practice, the number of columns of HRV seems to be 1.5x the one // of VISIR bands, but derive it from the paquet header poDS->nRasterXSize = nRawHRVColumns; }
CPLDebug("MSGN", "Deal with split %s", poDS->m_bHRVDealWithSplit == HRV_RSS ? "RSS" :CPLDebug("MSGN", "Deal with split %s", poDS->m_eHRVDealWithSplit == HRV_RSS ? "RSS" :else{ // In practice, the number of columns of HRV seems to be 1.5x the one // of VISIR bands, but derive it from the paquet header poDS->nRasterXSize = nRawHRVColumns; }
CPLDebug("MSGN", "Deal with split %s", poDS->m_bHRVDealWithSplit == HRV_RSS ? "RSS" :(poDS->m_bHRVDealWithSplit ? "whole" : "no"));(poDS->m_eHRVDealWithSplit ? "whole" : "no"));-- Reply to this email directly or view it on GitHub: https://github.com/OSGeo/gdal/pull/5866#pullrequestreview-995988632 You are receiving this because you authored the thread.
Message ID: @.***>
-- g8sqh @.***>
I've tested this, and I observe that within a RSS dataset, the georeferencing of the HRV and non-HRV band is consistent (when displaying in QGIS for example). But it is not consistent with MSGN datasets of whole-Earth images. There's a vertical shift of about a Earth radius. Some tweaking of the origin_y would be needed, or we'd probably need to declare the dataset height as being the whole Earth and tweak IReadBlock() to apply a vertical shift
Ok, the vertical shift is interesting (and I hadn't spotted it). At present I can't get qgis to convert from the geostationary projection (EPSG 9001) to lat/long - when I set the project CRS to (say) WGS-84 it won't show the image. (And I need to raise a QGIS question - I can't see how to pass HRV:
If qgis is able to read and convert EPSG 9001, then we need to keep the image location within the raster aligned with the projection data in the .nat file. gdalinfo on HRV:MSG4... looks strange - the origin is 5565, -5565km and the bottom right corner is (3x5565, -3x5565) I think this may be your one earth radius offset. I need to reread the nat file specification to understand this.
I can't see how to pass HRV: in
if you manually specify HRV:full_path_name in Add Raster layer, that works. But you can also use gdal_translate HRV:.... out.tif (or to a VRT file)
/home/travis/build/OSGeo/gdal/frmts/msgn/msgndataset.cpp: In static member function ‘static GDALDataset* MSGNDataset::Open(GDALOpenInfo*)’:
/home/travis/build/OSGeo/gdal/frmts/msgn/msgndataset.cpp:597:11: error: declaration of ‘i’ shadows a previous local [-Werror=shadow]
int i;
^
/home/travis/build/OSGeo/gdal/frmts/msgn/msgndataset.cpp:539:18: note: shadowed declaration is here
unsigned int i;
^
I believe this now works for MSG4, RSS from MSG3 and IODC from MSG2 - all for normal and HRV channels. I've hand-tweaked small constants to get decent registration to country boundaries - but no guarantee on lat/long from this
In qgis this will show the circular originals (with no CRS set), and also reproject onto a normal CRS (e.g. EPSG 3857)
I'll try again - I probably don't have the git skills to get this to work - may just keep as a patch set locally