Skip to content

Commit

Permalink
Merge pull request OSGeo#8435 from rouault/netCDF_EMIT
Browse files Browse the repository at this point in the history 
netCDF: add support for EMIT products and their geolocation and geometry lookup table
  • Loading branch information
@rouault
rouault authored Sep 22, 2023
2 parents 35ab41d + a7ec6b4 commit b3ee679
Show file tree
Hide file tree
Showing 12 changed files with 981 additions and 55 deletions.
105 changes: 76 additions & 29 deletions apps/gdalmdimtranslate_lib.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -602,7 +602,7 @@ static bool ParseArraySpec(const std::string &arraySpec, std::string &srcName,
std::string &dstName, int &band,
std::vector<int> &anTransposedAxis,
std::string &viewExpr,
GDALExtendedDataType &outputType)
GDALExtendedDataType &outputType, bool &bResampled)
{
if (!STARTS_WITH(arraySpec.c_str(), "name=") &&
!STARTS_WITH(arraySpec.c_str(), "band="))
Expand Down Expand Up @@ -695,6 +695,10 @@ static bool ParseArraySpec(const std::string &arraySpec, std::string &srcName,
outputType = GDALExtendedDataType::Create(eDT);
}
}
else if (STARTS_WITH(token.c_str(), "resample="))
{
bResampled = CPLTestBool(token.c_str() + strlen("resample="));
}
else
{
CPLError(CE_Failure, CPLE_AppDefined,
Expand All @@ -710,7 +714,9 @@ static bool ParseArraySpec(const std::string &arraySpec, std::string &srcName,
/************************************************************************/

static bool TranslateArray(
DimensionRemapper &oDimRemapper, const std::string &arraySpec,
DimensionRemapper &oDimRemapper,
const std::shared_ptr<GDALMDArray> &poSrcArrayIn,
const std::string &arraySpec,
const std::shared_ptr<GDALGroup> &poSrcRootGroup,
const std::shared_ptr<GDALGroup> &poSrcGroup,
const std::shared_ptr<GDALGroup> &poDstRootGroup,
Expand All @@ -724,14 +730,16 @@ static bool TranslateArray(
int band = -1;
std::vector<int> anTransposedAxis;
std::string viewExpr;
bool bResampled = false;
GDALExtendedDataType outputType(GDALExtendedDataType::Create(GDT_Unknown));
if (!ParseArraySpec(arraySpec, srcArrayName, dstArrayName, band,
anTransposedAxis, viewExpr, outputType))
anTransposedAxis, viewExpr, outputType, bResampled))
{
return false;
}

std::shared_ptr<GDALMDArray> srcArray;
bool bSrcArrayAccessibleThroughSrcGroup = true;
if (poSrcRootGroup && poSrcGroup)
{
if (!srcArrayName.empty() && srcArrayName[0] == '/')
Expand All @@ -740,9 +748,17 @@ static bool TranslateArray(
srcArray = poSrcGroup->OpenMDArray(srcArrayName);
if (!srcArray)
{
CPLError(CE_Failure, CPLE_AppDefined, "Cannot find array %s",
srcArrayName.c_str());
return false;
if (poSrcArrayIn && poSrcArrayIn->GetFullName() == arraySpec)
{
bSrcArrayAccessibleThroughSrcGroup = false;
srcArray = poSrcArrayIn;
}
else
{
CPLError(CE_Failure, CPLE_AppDefined, "Cannot find array %s",
srcArrayName.c_str());
return false;
}
}
}
else
Expand All @@ -754,6 +770,17 @@ static bool TranslateArray(
}

auto tmpArray = srcArray;

if (bResampled)
{
tmpArray =
tmpArray->GetResampled(std::vector<std::shared_ptr<GDALDimension>>(
tmpArray->GetDimensionCount()),
GRIORA_NearestNeighbour, nullptr, nullptr);
if (!tmpArray)
return false;
}

if (!anTransposedAxis.empty())
{
tmpArray = tmpArray->Transpose(anTransposedAxis);
Expand Down Expand Up @@ -1013,7 +1040,7 @@ static bool TranslateArray(
{
if (poSrcRootGroup)
{
if (!TranslateArray(oDimRemapper, indexingVarSpec,
if (!TranslateArray(oDimRemapper, srcIndexVar, indexingVarSpec,
poSrcRootGroup, poSrcGroup, poDstRootGroup,
poDstGroup, poSrcDS, mapSrcToDstDims,
mapDstDimFullNames, psOptions))
Expand Down Expand Up @@ -1072,6 +1099,8 @@ static bool TranslateArray(
GUInt64 nCurCost = 0;
dstArray->CopyFromAllExceptValues(srcArray.get(), false, nCurCost, 0,
nullptr, nullptr);
if (bResampled)
dstArray->SetSpatialRef(tmpArray->GetSpatialRef().get());

if (idxSliceSpec >= 0)
{
Expand Down Expand Up @@ -1139,23 +1168,39 @@ static bool TranslateArray(
}
}

const auto dimCount(tmpArray->GetDimensionCount());
std::vector<GUInt64> anSrcOffset(dimCount);
std::vector<GUInt64> anCount(dimCount);
for (size_t i = 0; i < dimCount; ++i)
double dfStart = 0.0;
double dfIncrement = 0.0;
if (!bSrcArrayAccessibleThroughSrcGroup &&
tmpArray->IsRegularlySpaced(dfStart, dfIncrement))
{
auto poSource = cpl::make_unique<VRTMDArraySourceRegularlySpaced>(
dfStart, dfIncrement);
dstArrayVRT->AddSource(std::move(poSource));
}
else
{
anCount[i] = tmpArrayDims[i]->GetSize();
const auto dimCount(tmpArray->GetDimensionCount());
std::vector<GUInt64> anSrcOffset(dimCount);
std::vector<GUInt64> anCount(dimCount);
for (size_t i = 0; i < dimCount; ++i)
{
anCount[i] = tmpArrayDims[i]->GetSize();
}
std::vector<GUInt64> anStep(dimCount, 1);
std::vector<GUInt64> anDstOffset(dimCount);
std::unique_ptr<VRTMDArraySourceFromArray> poSource(
new VRTMDArraySourceFromArray(
dstArrayVRT.get(), false, false, poSrcDS->GetDescription(),
band < 0 ? srcArray->GetFullName() : std::string(),
band >= 1 ? CPLSPrintf("%d", band) : std::string(),
std::move(anTransposedAxis),
bResampled ? (viewExpr.empty() ? std::string("resample=true")
: "resample=true," + viewExpr)
: viewExpr,
std::move(anSrcOffset), std::move(anCount), std::move(anStep),
std::move(anDstOffset)));
dstArrayVRT->AddSource(std::move(poSource));
}
std::vector<GUInt64> anStep(dimCount, 1);
std::vector<GUInt64> anDstOffset(dimCount);
std::unique_ptr<VRTMDArraySourceFromArray> poSource(
new VRTMDArraySourceFromArray(
dstArrayVRT.get(), false, false, poSrcDS->GetDescription(),
band < 0 ? srcArray->GetFullName() : std::string(),
band >= 1 ? CPLSPrintf("%d", band) : std::string(),
std::move(anTransposedAxis), viewExpr, std::move(anSrcOffset),
std::move(anCount), std::move(anStep), std::move(anDstOffset)));
dstArrayVRT->AddSource(std::move(poSource));

return true;
}
Expand Down Expand Up @@ -1245,8 +1290,8 @@ static bool CopyGroup(
auto arrayNames = poSrcGroup->GetMDArrayNames();
for (const auto &name : arrayNames)
{
if (!TranslateArray(oDimRemapper, name, poSrcRootGroup, poSrcGroup,
poDstRootGroup, poDstGroup, poSrcDS,
if (!TranslateArray(oDimRemapper, nullptr, name, poSrcRootGroup,
poSrcGroup, poDstRootGroup, poDstGroup, poSrcDS,
mapSrcToDstDims, mapDstDimFullNames, psOptions))
{
return false;
Expand Down Expand Up @@ -1431,10 +1476,10 @@ static bool TranslateInternal(std::shared_ptr<GDALGroup> &poDstRootGroup,
{
for (const auto &arraySpec : psOptions->aosArraySpec)
{
if (!TranslateArray(oDimRemapper, arraySpec, poSrcRootGroup,
poSrcRootGroup, poDstRootGroup, poDstRootGroup,
poSrcDS, mapSrcToDstDims, mapDstDimFullNames,
psOptions))
if (!TranslateArray(oDimRemapper, nullptr, arraySpec,
poSrcRootGroup, poSrcRootGroup, poDstRootGroup,
poDstRootGroup, poSrcDS, mapSrcToDstDims,
mapDstDimFullNames, psOptions))
{
return false;
}
Expand Down Expand Up @@ -1482,8 +1527,10 @@ CopyToNonMultiDimensionalDriver(GDALDriver *poDriver, const char *pszDest,
std::string viewExpr;
GDALExtendedDataType outputType(
GDALExtendedDataType::Create(GDT_Unknown));
bool bResampled = false;
ParseArraySpec(psOptions->aosArraySpec[0], srcArrayName, dstArrayName,
band, anTransposedAxis, viewExpr, outputType);
band, anTransposedAxis, viewExpr, outputType,
bResampled);
srcArray = poRG->OpenMDArray(dstArrayName);
}
else
Expand Down
Binary file added autotest/gdrivers/data/netcdf/fake_EMIT.nc
View file
Binary file not shown.
33 changes: 33 additions & 0 deletions autotest/gdrivers/netcdf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6418,3 +6418,36 @@ def test_netcdf_proj4string_geospatial_bounds_crs():
(-5400000.0, 75000.0, 0.0, 5400000.0, 0.0, -75000.0)
)
assert ds.GetSpatialRef().GetAuthorityCode(None) == "6931"


###############################################################################
# test opening a NASA EMIT dataset and check we get correct dimension mapping
# and a geolocation array


def test_netcdf_NASA_EMIT():

if not gdaltest.netcdf_drv_has_nc4:
pytest.skip("Requires NC4 support")

# Original dataset is https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20220903T163129_2224611_012/EMIT_L2A_RFL_001_20220903T163129_2224611_012.nc
ds = gdal.Open('NETCDF:"data/netcdf/fake_EMIT.nc":reflectance')
assert ds.RasterXSize == 2
assert ds.RasterYSize == 2
assert ds.RasterCount == 2
assert ds.GetRasterBand(1).ReadRaster() == struct.pack("f" * 4, 30, 40, 10, 20)
assert ds.GetRasterBand(2).ReadRaster() == struct.pack("f" * 4, -30, -40, -10, -20)

md = ds.GetMetadata("GEOLOCATION")
assert md == {
"GEOREFERENCING_CONVENTION": "PIXEL_CENTER",
"LINE_OFFSET": "0",
"LINE_STEP": "1",
"PIXEL_OFFSET": "0",
"PIXEL_STEP": "1",
"SRS": 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]',
"X_BAND": "1",
"X_DATASET": 'NETCDF:"data/netcdf/fake_EMIT.nc":/location/lon',
"Y_BAND": "1",
"Y_DATASET": 'NETCDF:"data/netcdf/fake_EMIT.nc":/location/lat',
}
Loading

0 comments on commit b3ee679

Please sign in to comment.