Skip to content

Commit

Permalink
Merge pull request #10215 from rouault/fix_10214
Browse files Browse the repository at this point in the history 
GRIB: make .idx reading compatible of /vsisubfile/
  • Loading branch information
@rouault
rouault committed Jun 18, 2024
2 parents 5818dea + 257fb20 commit da83066
Show file tree
Hide file tree
Showing 2 changed files with 105 additions and 30 deletions.
42 changes: 42 additions & 0 deletions autotest/gdrivers/grib.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2252,6 +2252,48 @@ def test_grib_grib2_sidecar():
) == ds_idx.GetRasterBand(i).GetMetadataItem(key)


def test_grib_grib2_sidecar_vsisubfile():

ds = gdal.Open("/vsisubfile/0_5359,data/grib/gfs.t06z.pgrb2.10p0.f010.grib2")
assert ds.RasterCount == 1
assert ds.GetRasterBand(1).GetDescription() == "REFD:1 hybrid level:10 hour fcst"

ds_ref = gdal.OpenEx(
"/vsisubfile/0_5359,data/grib/gfs.t06z.pgrb2.10p0.f010.grib2",
open_options=["USE_IDX=NO"],
)
assert ds_ref.RasterCount == 1
assert ds_ref.GetRasterBand(1).GetDescription() == '1[-] HYBL="Hybrid level"'
assert ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()

size = 16077 - 5359
ds = gdal.Open(f"/vsisubfile/5359_{size},data/grib/gfs.t06z.pgrb2.10p0.f010.grib2")
assert ds.RasterCount == 2
assert ds.GetRasterBand(1).GetDescription() == "REFD:2 hybrid level:10 hour fcst"
assert ds.GetRasterBand(2).GetDescription() == "REFC:entire atmosphere:10 hour fcst"

ds_ref = gdal.OpenEx(
f"/vsisubfile/5359_{size},data/grib/gfs.t06z.pgrb2.10p0.f010.grib2",
open_options=["USE_IDX=NO"],
)
assert ds_ref.RasterCount == 2
assert ds_ref.GetRasterBand(1).GetDescription() == '2[-] HYBL="Hybrid level"'
assert ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
assert ds.GetRasterBand(2).Checksum() == ds_ref.GetRasterBand(2).Checksum()

ds = gdal.Open("/vsisubfile/16077_-1,data/grib/gfs.t06z.pgrb2.10p0.f010.grib2")
assert ds.RasterCount == 3
assert ds.GetRasterBand(1).GetDescription() == "VIS:surface:10 hour fcst"
assert (
ds.GetRasterBand(2).GetDescription()
== "UGRD:planetary boundary layer:10 hour fcst"
)
assert (
ds.GetRasterBand(3).GetDescription()
== "VGRD:planetary boundary layer:10 hour fcst"
)


# Test reading a (broken) mix of GRIBv2/GRIBv1 bands


Expand Down
93 changes: 63 additions & 30 deletions frmts/grib/gribdataset.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -922,11 +922,16 @@ char **GRIBRasterBand::GetMetadata(const char *pszDomain)
const char *GRIBRasterBand::GetMetadataItem(const char *pszName,
const char *pszDomain)
{
FindMetaData();
if (m_nGribVersion == 2 &&
CPLTestBool(CPLGetConfigOption("GRIB_PDS_ALL_BANDS", "ON")))
if (!((!pszDomain || pszDomain[0] == 0) &&
(EQUAL(pszName, "STATISTICS_MINIMUM") ||
EQUAL(pszName, "STATISTICS_MAXIMUM"))))
{
FindPDSTemplateGRIB2();
FindMetaData();
if (m_nGribVersion == 2 &&
CPLTestBool(CPLGetConfigOption("GRIB_PDS_ALL_BANDS", "ON")))
{
FindPDSTemplateGRIB2();
}
}
return GDALPamRasterBand::GetMetadataItem(pszName, pszDomain);
}
Expand Down Expand Up @@ -1156,34 +1161,34 @@ class InventoryWrapperGrib : public gdal::grib::InventoryWrapper
class InventoryWrapperSidecar : public gdal::grib::InventoryWrapper
{
public:
explicit InventoryWrapperSidecar(VSILFILE *fp)
explicit InventoryWrapperSidecar(VSILFILE *fp, uint64_t nStartOffset,
int64_t nSize)
: gdal::grib::InventoryWrapper()
{
result_ = -1;
VSIFSeekL(fp, 0, SEEK_END);
size_t length = static_cast<size_t>(VSIFTellL(fp));
if (length > 4 * 1024 * 1024)
return;
std::string psSidecar;
psSidecar.resize(length);
std::string osSidecar;
osSidecar.resize(length);
VSIFSeekL(fp, 0, SEEK_SET);
if (VSIFReadL(&psSidecar[0], length, 1, fp) != 1)
if (VSIFReadL(&osSidecar[0], length, 1, fp) != 1)
return;

CPLStringList aosMsgs(
CSLTokenizeString2(psSidecar.c_str(), "\n",
const CPLStringList aosMsgs(
CSLTokenizeString2(osSidecar.c_str(), "\n",
CSLT_PRESERVEQUOTES | CSLT_STRIPLEADSPACES));
inv_len_ = aosMsgs.size();
inv_ = static_cast<inventoryType *>(
CPLMalloc(inv_len_ * sizeof(inventoryType)));
CPLCalloc(aosMsgs.size(), sizeof(inventoryType)));

for (size_t i = 0; i < inv_len_; ++i)
for (const char *pszMsg : aosMsgs)
{
// We are parsing
// "msgNum[.subgNum]:start:dontcare:name1:name2:name3" For NOMADS:
// "msgNum[.subgNum]:start:reftime:var:level:time"
CPLStringList aosTokens(CSLTokenizeString2(
aosMsgs[i], ":", CSLT_PRESERVEQUOTES | CSLT_ALLOWEMPTYTOKENS));
const CPLStringList aosTokens(CSLTokenizeString2(
pszMsg, ":", CSLT_PRESERVEQUOTES | CSLT_ALLOWEMPTYTOKENS));
CPLStringList aosNum;

if (aosTokens.size() < 6)
Expand All @@ -1200,7 +1205,7 @@ class InventoryWrapperSidecar : public gdal::grib::InventoryWrapper
goto err_sidecar;

if (aosNum.size() < 2)
inv_[i].subgNum = 0;
inv_[inv_len_].subgNum = 0;
else
{
auto subgNum = strtol(aosNum[1], &endptr, 10);
Expand All @@ -1211,30 +1216,37 @@ class InventoryWrapperSidecar : public gdal::grib::InventoryWrapper
// .idx file use a 1-based indexing, whereas DEGRIB uses a
// 0-based one
subgNum--;
inv_[i].subgNum = static_cast<unsigned short>(subgNum);
inv_[inv_len_].subgNum = static_cast<unsigned short>(subgNum);
}

inv_[i].start = strtoll(aosTokens[1], &endptr, 10);
inv_[inv_len_].start = strtoll(aosTokens[1], &endptr, 10);
if (*endptr != 0)
goto err_sidecar;

inv_[i].unitName = nullptr;
inv_[i].comment = nullptr;
inv_[i].element = nullptr;
inv_[i].shortFstLevel = nullptr;
if (inv_[inv_len_].start < nStartOffset)
continue;
if (nSize > 0 && inv_[inv_len_].start >= nStartOffset + nSize)
break;

inv_[inv_len_].start -= nStartOffset;

inv_[inv_len_].unitName = nullptr;
inv_[inv_len_].comment = nullptr;
inv_[inv_len_].element = nullptr;
inv_[inv_len_].shortFstLevel = nullptr;
// This is going into the description field ->
// the only one available before loading the metadata
inv_[i].longFstLevel = VSIStrdup(CPLSPrintf(
inv_[inv_len_].longFstLevel = VSIStrdup(CPLSPrintf(
"%s:%s:%s", aosTokens[3], aosTokens[4], aosTokens[5]));
++inv_len_;

continue;

err_sidecar:
CPLDebug("GRIB",
"Failed parsing sidecar entry '%s', "
"falling back to constructing an inventory",
aosMsgs[i]);
inv_len_ = static_cast<unsigned>(i);
pszMsg);
return;
}

Expand Down Expand Up @@ -1309,22 +1321,43 @@ GRIBDataset::Inventory(VSILFILE *fp, GDALOpenInfo *poOpenInfo)
std::unique_ptr<gdal::grib::InventoryWrapper> pInventories;

VSIFSeekL(fp, 0, SEEK_SET);
CPLString sSideCarFilename = CPLString(poOpenInfo->pszFilename) + ".idx";
std::string osSideCarFilename(poOpenInfo->pszFilename);
uint64_t nStartOffset = 0;
int64_t nSize = -1;
if (STARTS_WITH(poOpenInfo->pszFilename, "/vsisubfile/"))
{
const char *pszPtr = poOpenInfo->pszFilename + strlen("/vsisubfile/");
const char *pszComma = strchr(pszPtr, ',');
if (pszComma)
{
const CPLStringList aosTokens(CSLTokenizeString2(
std::string(pszPtr, pszComma - pszPtr).c_str(), "_", 0));
if (aosTokens.size() == 2)
{
nStartOffset = std::strtoull(aosTokens[0], nullptr, 10);
nSize = std::strtoll(aosTokens[1], nullptr, 10);
osSideCarFilename = pszComma + 1;
}
}
}
osSideCarFilename += ".idx";
VSILFILE *fpSideCar = nullptr;
if (CPLTestBool(CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
"USE_IDX", "YES")) &&
((fpSideCar = VSIFOpenL(sSideCarFilename, "rb")) != nullptr))
((fpSideCar = VSIFOpenL(osSideCarFilename.c_str(), "rb")) != nullptr))
{
CPLDebug("GRIB", "Reading inventories from sidecar file %s",
sSideCarFilename.c_str());
osSideCarFilename.c_str());
// Contains an GRIB2 message inventory of the file.
pInventories = std::make_unique<InventoryWrapperSidecar>(fpSideCar);
pInventories = std::make_unique<InventoryWrapperSidecar>(
fpSideCar, nStartOffset, nSize);
if (pInventories->result() <= 0 || pInventories->length() == 0)
pInventories = nullptr;
VSIFCloseL(fpSideCar);
}
else
CPLDebug("GRIB", "Failed opening sidecar %s", sSideCarFilename.c_str());
CPLDebug("GRIB", "Failed opening sidecar %s",
osSideCarFilename.c_str());

if (pInventories == nullptr)
{
Expand Down

0 comments on commit da83066

Please sign in to comment.