From cd9638b9c2dcc0224feb5b839d14508e9dc8835d Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 19 Sep 2023 23:03:39 +0200 Subject: [PATCH 1/6] netCDF (traditional): add support for EMIT band data ordering and geolocation array (not using glt_x/glt_y) --- autotest/gdrivers/data/netcdf/fake_EMIT.nc | Bin 0 -> 12164 bytes autotest/gdrivers/netcdf.py | 33 ++++++ frmts/netcdf/netcdfdataset.cpp | 127 ++++++++++++++++++++- frmts/netcdf/netcdfdataset.h | 2 + 4 files changed, 160 insertions(+), 2 deletions(-) create mode 100644 autotest/gdrivers/data/netcdf/fake_EMIT.nc diff --git a/autotest/gdrivers/data/netcdf/fake_EMIT.nc b/autotest/gdrivers/data/netcdf/fake_EMIT.nc new file mode 100644 index 0000000000000000000000000000000000000000..4888352fc5622e0db19c05aca55d2c8a27e03fe1 GIT binary patch literal 12164 zcmeGiS!^3cbT&TfrkgZ{LJI`gx&;aiCb6Bv&_aXlB&JRrVy7Hc`iIE z;G^P@kT|{(hk$~F%11vC;scd}xFnzgQ4Vo^a3e&33W0Wp| z#?FNw4!|b_810Xx#JnAKP=#_xHo2RzrZ|amamKpCwDbDN&OaU4!q^K(-db9*RU8@% zEv6FER^8_#l|tmqk+V!&u!~U}P3*|&cHbJLXUx44|0_w{&}%6p@U-f(h@ z_m0e=U!3=YNmqcB@<-Z0ENvLidZEv;3^yF#sP zNIpCeOKtIHj4>bwN?&VRXRtls5Bb~MLZM)w+8RP_em`HstufeS9C%?*W}x#xVkp%I zOA~0z9FzEHcQRFN#X+%R^;X8v5*mQJE7_P~fwGx}Q7b==%xwZQ5bB~~l_`qgva}XP z5V6F6QZLj(*E0KeuB;846-8)~9D#I+pkp}>ux!<(#ralv*or%0T$cx1>gwwF1ZXzk zIj+LhF!qlCMFVxuXaj01k}9xbG|+IoJDG$!>JHjJtvv$L8%V|y1LvwHrw0ssBr_U}WXH1^m}Qy_09PvO4!fL!!yue_;jrf%9NnRi7a>nbPH?_? zEF9qhsY=hvtqE-7I#K;rVnIb`DQl*lx)=S0=sPMG-b+1P^uXK~<73I8)sN(@0=R5& zv$spC{RbZ1%6-qag{uLmfA60A@SZa(|5;vwr(|P7TKlE6j+g)Qt58sky)FQoUH;A= z`cfz}_eDrvde3=s)(H4X%iJU_R!pgjP=xi+;L;au|7-v1g)VwrYU;*yi>iyBUd2Ut zg3h(l^i|0jGJV%PyH_EHXQUjC7romjCO<;bnbq_k-dBHr$@I6#g;z?}eIX8}IO!F& z$jaM0tn{AY#bS~Vo?h|6nUZ|q6x}|`)Ot6u^Zy|EC#0N?m;6!b^iNFkiNHtuN|GOL z6q3A20DQ0VLXxL5`=aR{eIfz*jY8E|P; z0xlIJ8w|L(ovd$u-N{2Rs*#uYOj@Xh%YZVsmk{>;n3T~@>5bL35N^PfA~c&cvvvVS z!_dqFKS{(8pEC*t86QH+FBj~aDoOrr{sqR^E@qO94&qGuAk{e5#`nVjBa}PO4J(838QGLYF zhDw`PCoqy$*iwjqSF%#%8;>P{4HrOXjEq;nkBec3P2vRgNC*WLq`Z#CR$fODyz}52 zzxRd{alSCwk3T~j&Vz+MbxD;g@x*AtQ_RK{6#5gM#v!98JkgU%6!9`=K&b=5(>PQo zNO<}7D9%gc%lBFlPxR8~Ov6iWSj5vfQcH_?dR<97nucG>xme-@IXTJhY4pn~=xKcZ zAyCFo8b^||jkEwl)L@MRH4fA`P~$+212qoRI8ft2jRQ3f)HqP%K#c=s9KciyC6U(4 z)BZKG9Z%d&J}W_)0)pjP2mDAaujAM8Fd3GD*e5QAtCbs?vp{rk#JvOcd@P3qeZekq*3i;N@14Dt%K?)^bX5giJZ%86l-c zZ;+XijV!o?1&=}6NhC1uzmiYP_6oq9JYh=vK#7xf|LqebnxvdNMXC-Wsd9=bNZLC= zF&3mBVbs~ukJER`d(OaVY~y-;?OZ8H%pj3`;&dlLTMJh+uI>t6iI--oCw@_%sh_wV zQP2#%V6S$pmZ#nGEz_?wm=em|GHsR~OXKt>lJ-8UWEz@FrkE`x+0gZql1Q(fNu-#> zq;4cq6VKm5w;X&#*y;WXJEc7DjF8)1xqomScXaXV$H$J9jxMk!?6+Qu?(Y9cg18@i z{iqRS2_6dI-J8lof%1LfLT5S*zEP#GB*rY$Zk+@@vM1zV*8_E?{8-^R>Aio zZh%gS@U-N};h!VrQDEg=BlTmr_XWKcX5>-nk0|nhjqyTvDydHQD5<{qyBF8<`0A|O zB`G8AlIpZes(;aZ-j95jZ&F0R5Ti(nlnBzTR@w#KW2O4m_pi5cS9B+qbiuV+ZU7k)J6E+dbZ1vmIx;Qk<7{o4@qb@~D= zJ2Io~8`^wrt-j#VrcHDhM@n?ku-x9Z$rxsh9A*B$P$)_utHHY;cotOXs literal 0 HcmV?d00001 diff --git a/autotest/gdrivers/netcdf.py b/autotest/gdrivers/netcdf.py index 19ef9314d2af..d8b9c414c781 100755 --- a/autotest/gdrivers/netcdf.py +++ b/autotest/gdrivers/netcdf.py @@ -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', + } diff --git a/frmts/netcdf/netcdfdataset.cpp b/frmts/netcdf/netcdfdataset.cpp index 9e6c1732ff50..91e6fa437ae7 100644 --- a/frmts/netcdf/netcdfdataset.cpp +++ b/frmts/netcdf/netcdfdataset.cpp @@ -5283,6 +5283,114 @@ bool netCDFDataset::ProcessNASAL2OceanGeoLocation(int nGroupId, int nVarId) "GEOLOCATION"); return true; } + +bool netCDFDataset::ProcessNASAEMITGeoLocation(int nGroupId, int nVarId) +{ + // Cf https://earth.jpl.nasa.gov/emit/data/data-portal/coverage-and-forecasts/ + + // Check for a structure like: + /* netcdf EMIT_L2A_RFL_001_20220903T163129_2224611_012 { + dimensions: + downtrack = 1280 ; + crosstrack = 1242 ; + bands = 285 ; + [...] + + variables: + float reflectance(downtrack, crosstrack, bands) ; + + group: location { + variables: + double lon(downtrack, crosstrack) ; + lon:_FillValue = -9999. ; + lon:long_name = "Longitude (WGS-84)" ; + lon:units = "degrees east" ; + double lat(downtrack, crosstrack) ; + lat:_FillValue = -9999. ; + lat:long_name = "Latitude (WGS-84)" ; + lat:units = "degrees north" ; + } // group location + + } + */ + + int nVarDims = 0; + NCDF_ERR(nc_inq_varndims(nGroupId, nVarId, &nVarDims)); + if (nVarDims != 3) + return false; + + int nLocationGrpId = 0; + if (nc_inq_grp_ncid(cdfid, "location", &nLocationGrpId) != NC_NOERR) + return false; + + std::array anVarDimIds; + NCDF_ERR(nc_inq_vardimid(nGroupId, nVarId, anVarDimIds.data())); + if (nYDimID != anVarDimIds[0] || nXDimID != anVarDimIds[1]) + return false; + + int nLongitudeId = 0; + int nLatitudeId = 0; + if (nc_inq_varid(nLocationGrpId, "lon", &nLongitudeId) != NC_NOERR || + nc_inq_varid(nLocationGrpId, "lat", &nLatitudeId) != NC_NOERR) + { + return false; + } + + int nDimsLongitude = 0; + NCDF_ERR(nc_inq_varndims(nLocationGrpId, nLongitudeId, &nDimsLongitude)); + int nDimsLatitude = 0; + NCDF_ERR(nc_inq_varndims(nLocationGrpId, nLatitudeId, &nDimsLatitude)); + if (!(nDimsLongitude == 2 && nDimsLatitude == 2)) + { + return false; + } + + std::array anDimLongitudeIds; + NCDF_ERR(nc_inq_vardimid(nLocationGrpId, nLongitudeId, + anDimLongitudeIds.data())); + std::array anDimLatitudeIds; + NCDF_ERR( + nc_inq_vardimid(nLocationGrpId, nLatitudeId, anDimLatitudeIds.data())); + if (anDimLongitudeIds != anDimLatitudeIds) + { + return false; + } + + if (anDimLongitudeIds[0] != anVarDimIds[0] || + anDimLongitudeIds[1] != anVarDimIds[1]) + { + return false; + } + + const char *pszGeolocXFullName = "/location/lon"; + const char *pszGeolocYFullName = "/location/lat"; + + CPLDebug("GDAL_netCDF", "using variables %s and %s for GEOLOCATION", + pszGeolocXFullName, pszGeolocYFullName); + + GDALPamDataset::SetMetadataItem("SRS", SRS_WKT_WGS84_LAT_LONG, + "GEOLOCATION"); + + CPLString osTMP; + osTMP.Printf("NETCDF:\"%s\":%s", osFilename.c_str(), pszGeolocXFullName); + + GDALPamDataset::SetMetadataItem("X_DATASET", osTMP, "GEOLOCATION"); + GDALPamDataset::SetMetadataItem("X_BAND", "1", "GEOLOCATION"); + osTMP.Printf("NETCDF:\"%s\":%s", osFilename.c_str(), pszGeolocYFullName); + + GDALPamDataset::SetMetadataItem("Y_DATASET", osTMP, "GEOLOCATION"); + GDALPamDataset::SetMetadataItem("Y_BAND", "1", "GEOLOCATION"); + + GDALPamDataset::SetMetadataItem("PIXEL_OFFSET", "0", "GEOLOCATION"); + GDALPamDataset::SetMetadataItem("PIXEL_STEP", "1", "GEOLOCATION"); + + GDALPamDataset::SetMetadataItem("LINE_OFFSET", "0", "GEOLOCATION"); + GDALPamDataset::SetMetadataItem("LINE_STEP", "1", "GEOLOCATION"); + + GDALPamDataset::SetMetadataItem("GEOREFERENCING_CONVENTION", "PIXEL_CENTER", + "GEOLOCATION"); + return true; +} #endif int netCDFDataset::ProcessCFGeolocation(int nGroupId, int nVarId, @@ -5423,6 +5531,9 @@ int netCDFDataset::ProcessCFGeolocation(int nGroupId, int nVarId, else { bAddGeoloc = ProcessNASAL2OceanGeoLocation(nGroupId, nVarId); + + if (!bAddGeoloc) + bAddGeoloc = ProcessNASAEMITGeoLocation(nGroupId, nVarId); } #endif @@ -9594,16 +9705,28 @@ GDALDataset *netCDFDataset::Open(GDALOpenInfo *poOpenInfo) } } + // For example for EMIT data (https://earth.jpl.nasa.gov/emit/data/data-portal/coverage-and-forecasts/), + // dimension order is downtrack, crosstrack, bands + bool bYXBandOrder = false; + if (nd == 3) + { + char szDimName[NC_MAX_NAME + 1] = {}; + status = nc_inq_dimname(cdfid, poDS->m_anDimIds[2], szDimName); + NCDF_ERR(status); + bYXBandOrder = + strcmp(szDimName, "bands") == 0 || strcmp(szDimName, "band") == 0; + } + // Get X dimensions information. size_t xdim; - poDS->nXDimID = poDS->m_anDimIds[nd - 1]; + poDS->nXDimID = poDS->m_anDimIds[bYXBandOrder ? 1 : nd - 1]; nc_inq_dimlen(cdfid, poDS->nXDimID, &xdim); // Get Y dimension information. size_t ydim; if (nd >= 2) { - poDS->nYDimID = poDS->m_anDimIds[nd - 2]; + poDS->nYDimID = poDS->m_anDimIds[bYXBandOrder ? 0 : nd - 2]; nc_inq_dimlen(cdfid, poDS->nYDimID, &ydim); } else diff --git a/frmts/netcdf/netcdfdataset.h b/frmts/netcdf/netcdfdataset.h index 3e988061fc2a..e2a4a81d02f5 100644 --- a/frmts/netcdf/netcdfdataset.h +++ b/frmts/netcdf/netcdfdataset.h @@ -907,6 +907,8 @@ class netCDFDataset final : public GDALPamDataset #ifdef NETCDF_HAS_NC4 bool ProcessNASAL2OceanGeoLocation(int nGroupId, int nVarId); + + bool ProcessNASAEMITGeoLocation(int nGroupId, int nVarId); #endif int ProcessCFGeolocation(int nGroupId, int nVarId, From 8ddeba0551dce5efa3963882d2946b5ed89144b5 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 20 Sep 2023 03:05:11 +0200 Subject: [PATCH 2/6] GDALMDArrayResampled::Create(): fix assertion if providing a input non-empty dimension for the non-X/Y dimensions --- gcore/gdalmultidim.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/gcore/gdalmultidim.cpp b/gcore/gdalmultidim.cpp index c85ab2789c7c..7b6337e8ab8d 100644 --- a/gcore/gdalmultidim.cpp +++ b/gcore/gdalmultidim.cpp @@ -7758,6 +7758,10 @@ std::shared_ptr GDALMDArrayResampled::Create( i); return nullptr; } + else + { + apoNewDims.emplace_back(aoParentDims[i]); + } anBlockSize.emplace_back(anParentBlockSize[i]); } From c60d9d663d3ef8fd6fca88923265979a2539f515 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 20 Sep 2023 01:59:03 +0200 Subject: [PATCH 3/6] netCDF multidimensional: add support for EMIT band data ordering and geolocation array (not using glt_x/glt_y) --- autotest/gdrivers/netcdf_multidim.py | 70 ++++++++++++++++++++++++++++ frmts/netcdf/netcdfmultidim.cpp | 22 +++++++++ gcore/gdalmultidim.cpp | 60 +++++++++++++++++------- 3 files changed, 135 insertions(+), 17 deletions(-) diff --git a/autotest/gdrivers/netcdf_multidim.py b/autotest/gdrivers/netcdf_multidim.py index 83eca0dd791e..ab0721328e0f 100755 --- a/autotest/gdrivers/netcdf_multidim.py +++ b/autotest/gdrivers/netcdf_multidim.py @@ -3576,3 +3576,73 @@ def reopen(): finally: gdal.Unlink(filename) gdal.Unlink(filename + ".aux.xml") + + +def test_netcdf_multidim_getresampled_with_geoloc_EMIT(): + + ds = gdal.OpenEx("data/netcdf/fake_EMIT.nc", gdal.OF_MULTIDIM_RASTER) + rg = ds.GetRootGroup() + + ar = rg.OpenMDArray("reflectance") + coordinate_vars = ar.GetCoordinateVariables() + assert len(coordinate_vars) == 2 + assert coordinate_vars[0].GetName() == "lon" + assert coordinate_vars[1].GetName() == "lat" + + resampled_ar = ar.GetResampled( + [None, None, ar.GetDimensions()[2]], gdal.GRIORA_NearestNeighbour, None + ) + assert resampled_ar is not None + dims = resampled_ar.GetDimensions() + assert dims[0].GetName() == "dimY" + assert dims[0].GetSize() == 3 + assert dims[1].GetName() == "dimX" + assert dims[1].GetSize() == 3 + assert dims[2].GetName() == "bands" + assert dims[2].GetSize() == 2 + + resampled_ar = ar.GetResampled( + [None] * ar.GetDimensionCount(), gdal.GRIORA_NearestNeighbour, None + ) + assert resampled_ar is not None + dims = resampled_ar.GetDimensions() + assert dims[0].GetName() == "dimY" + assert dims[0].GetSize() == 3 + assert dims[1].GetName() == "dimX" + assert dims[1].GetSize() == 3 + assert dims[2].GetName() == "bands" + assert dims[2].GetSize() == 2 + + resampled_ar_transposed = resampled_ar.Transpose([2, 0, 1]) + dims = resampled_ar_transposed.GetDimensions() + assert dims[0].GetName() == "bands" + assert dims[0].GetSize() == 2 + assert dims[1].GetName() == "dimY" + assert dims[1].GetSize() == 3 + assert dims[2].GetName() == "dimX" + assert dims[2].GetSize() == 3 + + # By default, the classic netCDF driver would use bottom-up reordering, + # which slightly modifies the output of the geolocation interpolation, + # and would not make it possible to compare exactly with the GetResampled() + # result + with gdaltest.config_option("GDAL_NETCDF_BOTTOMUP", "NO"): + warped_ds = gdal.Warp( + "", 'NETCDF:"data/netcdf/fake_EMIT.nc":reflectance', format="MEM" + ) + assert warped_ds.ReadRaster() == resampled_ar_transposed.Read() + xoff = 1 + yoff = 2 + xsize = 2 + ysize = 1 + band_count = 2 + assert warped_ds.ReadRaster( + xoff, yoff, xsize, ysize + ) == resampled_ar_transposed.Read( + array_start_idx=[0, yoff, xoff], count=[band_count, ysize, xsize] + ) + assert warped_ds.GetRasterBand(2).ReadRaster( + xoff, yoff, xsize, ysize + ) == resampled_ar_transposed.Read( + array_start_idx=[1, yoff, xoff], count=[1, ysize, xsize] + ) diff --git a/frmts/netcdf/netcdfmultidim.cpp b/frmts/netcdf/netcdfmultidim.cpp index 574cc201cf0a..f2c0e10c2334 100644 --- a/frmts/netcdf/netcdfmultidim.cpp +++ b/frmts/netcdf/netcdfmultidim.cpp @@ -4040,6 +4040,28 @@ netCDFVariable::GetCoordinateVariables() const } } + // Special case for NASA EMIT datasets + auto apoDims = GetDimensions(); + if (apoDims.size() == 3 && apoDims[0]->GetName() == "downtrack" && + apoDims[1]->GetName() == "crosstrack" && + apoDims[2]->GetName() == "bands") + { + auto poRootGroup = netCDFGroup::Create(m_poShared, nullptr, m_gid); + if (poRootGroup) + { + auto poLocationGroup = poRootGroup->OpenGroup("location"); + if (poLocationGroup) + { + auto poLon = poLocationGroup->OpenMDArray("lon"); + auto poLat = poLocationGroup->OpenMDArray("lat"); + if (poLon && poLat) + { + return {poLon, poLat}; + } + } + } + } + return ret; } diff --git a/gcore/gdalmultidim.cpp b/gcore/gdalmultidim.cpp index 7b6337e8ab8d..ec37cd92458e 100644 --- a/gcore/gdalmultidim.cpp +++ b/gcore/gdalmultidim.cpp @@ -7366,8 +7366,8 @@ class GDALMDArrayResampledDataset final : public GDALPamDataset friend class GDALMDArrayResampledDatasetRasterBand; std::shared_ptr m_poArray; - size_t m_iXDim; - size_t m_iYDim; + const size_t m_iXDim; + const size_t m_iYDim; double m_adfGeoTransform[6]{0, 1, 0, 0, 0, 1}; bool m_bHasGT = false; mutable std::shared_ptr m_poSRS{}; @@ -7586,7 +7586,7 @@ class GDALMDArrayResampled final : public GDALPamMDArray void *pDstBuffer) const override; public: - static std::shared_ptr + static std::shared_ptr Create(const std::shared_ptr &poParent, const std::vector> &apoNewDims, GDALRIOResampleAlg resampleAlg, @@ -7669,7 +7669,7 @@ class GDALMDArrayResampled final : public GDALPamMDArray /* GDALMDArrayResampled::Create() */ /************************************************************************/ -std::shared_ptr GDALMDArrayResampled::Create( +std::shared_ptr GDALMDArrayResampled::Create( const std::shared_ptr &poParent, const std::vector> &apoNewDimsIn, GDALRIOResampleAlg resampleAlg, const OGRSpatialReference *poTargetSRS, @@ -7743,8 +7743,22 @@ std::shared_ptr GDALMDArrayResampled::Create( anBlockSize.reserve(apoNewDimsIn.size()); const auto &anParentBlockSize = poParent->GetBlockSize(); - for (unsigned i = 0; i + 2 < apoNewDimsIn.size(); ++i) + auto apoParentDims = poParent->GetDimensions(); + // Special case for NASA EMIT datasets + const bool bYXBandOrder = (apoParentDims.size() == 3 && + apoParentDims[0]->GetName() == "downtrack" && + apoParentDims[1]->GetName() == "crosstrack" && + apoParentDims[2]->GetName() == "bands"); + + const size_t iYDimParent = + bYXBandOrder ? 0 : poParent->GetDimensionCount() - 2; + const size_t iXDimParent = + bYXBandOrder ? 1 : poParent->GetDimensionCount() - 1; + + for (unsigned i = 0; i < apoNewDimsIn.size(); ++i) { + if (i == iYDimParent || i == iXDimParent) + continue; if (apoNewDimsIn[i] == nullptr) { apoNewDims.emplace_back(aoParentDims[i]); @@ -7765,15 +7779,13 @@ std::shared_ptr GDALMDArrayResampled::Create( anBlockSize.emplace_back(anParentBlockSize[i]); } - const size_t iYDim = poParent->GetDimensionCount() - 2; - const size_t iXDim = poParent->GetDimensionCount() - 1; std::unique_ptr poParentDS( - new GDALMDArrayResampledDataset(poParent, iXDim, iYDim)); + new GDALMDArrayResampledDataset(poParent, iXDimParent, iYDimParent)); double dfXStart = 0.0; double dfXSpacing = 0.0; bool gotXSpacing = false; - auto poNewDimX = apoNewDimsIn[iXDim]; + auto poNewDimX = apoNewDimsIn[iXDimParent]; if (poNewDimX) { if (poNewDimX->GetSize() > static_cast(INT_MAX)) @@ -7800,7 +7812,7 @@ std::shared_ptr GDALMDArrayResampled::Create( double dfYStart = 0.0; double dfYSpacing = 0.0; - auto poNewDimY = apoNewDimsIn[iYDim]; + auto poNewDimY = apoNewDimsIn[iYDimParent]; bool gotYSpacing = false; if (poNewDimY) { @@ -7884,8 +7896,8 @@ std::shared_ptr GDALMDArrayResampled::Create( const auto &longDims = poLongVar->GetDimensions(); const auto latDimCount = poLatVar->GetDimensionCount(); const auto &latDims = poLatVar->GetDimensions(); - const auto xDimSize = aoParentDims[iXDim]->GetSize(); - const auto yDimSize = aoParentDims[iYDim]->GetSize(); + const auto xDimSize = aoParentDims[iXDimParent]->GetSize(); + const auto yDimSize = aoParentDims[iYDimParent]->GetSize(); if (longDimCount == 1 && longDims[0]->GetSize() == xDimSize && latDimCount == 1 && latDims[0]->GetSize() == yDimSize) { @@ -8074,6 +8086,14 @@ std::shared_ptr GDALMDArrayResampled::Create( newAr->m_poVarY = varY; newAr->m_poReprojectedDS = std::move(poReprojectedDS); newAr->m_poParentDS = std::move(poParentDS); + + // If the input array is y,x,band ordered, the above newAr is + // actually band,y,x ordered as it is more convenient for + // GDALMDArrayResampled::IRead() implementation. But transpose that + // array to the order of the input array + if (bYXBandOrder) + return newAr->Transpose(std::vector{1, 2, 0}); + return newAr; } @@ -8112,6 +8132,8 @@ bool GDALMDArrayResampled::IRead(const GUInt64 *arrayStartIdx, // Use an array to avoid a false positive warning from CLang Static // Analyzer about flushCaches being never read bool flushCaches[] = {false}; + const bool bYXBandOrder = + m_poParentDS->m_iYDim == 0 && m_poParentDS->m_iXDim == 1; lbl_next_depth: if (dimIdx == iDimY) @@ -8135,11 +8157,13 @@ bool GDALMDArrayResampled::IRead(const GUInt64 *arrayStartIdx, else { stack[dimIdx].nIters = count[dimIdx]; - if (m_poParentDS->m_anOffset[dimIdx] != arrayStartIdx[dimIdx]) + if (m_poParentDS->m_anOffset[bYXBandOrder ? 2 : dimIdx] != + arrayStartIdx[dimIdx]) { flushCaches[0] = true; } - m_poParentDS->m_anOffset[dimIdx] = arrayStartIdx[dimIdx]; + m_poParentDS->m_anOffset[bYXBandOrder ? 2 : dimIdx] = + arrayStartIdx[dimIdx]; while (true) { dimIdx++; @@ -8150,7 +8174,7 @@ bool GDALMDArrayResampled::IRead(const GUInt64 *arrayStartIdx, if ((--stack[dimIdx].nIters) == 0) break; flushCaches[0] = true; - ++m_poParentDS->m_anOffset[dimIdx]; + ++m_poParentDS->m_anOffset[bYXBandOrder ? 2 : dimIdx]; stack[dimIdx].dst_ptr += stack[dimIdx].dst_inc_offset; } } @@ -8168,7 +8192,8 @@ bool GDALMDArrayResampled::IRead(const GUInt64 *arrayStartIdx, * * This is the same as the C function GDALMDArrayGetResampled(). * - * Currently this method can only resample along the last 2 dimensions. + * Currently this method can only resample along the last 2 dimensions, unless + * orthorectifying a NASA EMIT dataset. * * @param apoNewDims New dimensions. Its size should be GetDimensionCount(). * apoNewDims[i] can be NULL to let the method automatically @@ -11647,7 +11672,8 @@ GDALMDArrayH GDALMDArrayGetMask(GDALMDArrayH hArray, CSLConstList papszOptions) * * This is the same as the C++ method GDALMDArray::GetResampled(). * - * Currently this method can only resample along the last 2 dimensions. + * Currently this method can only resample along the last 2 dimensions, unless + * orthorectifying a NASA EMIT dataset. * * The returned object should be released with GDALMDArrayRelease(). * From a7ac4e035b311b14c9e1ff36e66a57832409549f Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 20 Sep 2023 04:45:48 +0200 Subject: [PATCH 4/6] netCDF multidimensional: use glt_x/glt_y for EMIT data --- autotest/gdrivers/netcdf_multidim.py | 98 +++++- frmts/netcdf/netcdfmultidim.cpp | 427 +++++++++++++++++++++++++++ gcore/gdal_priv.h | 2 +- gcore/gdalmultidim.cpp | 3 + 4 files changed, 527 insertions(+), 3 deletions(-) diff --git a/autotest/gdrivers/netcdf_multidim.py b/autotest/gdrivers/netcdf_multidim.py index ab0721328e0f..06ac484ca6a0 100755 --- a/autotest/gdrivers/netcdf_multidim.py +++ b/autotest/gdrivers/netcdf_multidim.py @@ -3590,7 +3590,10 @@ def test_netcdf_multidim_getresampled_with_geoloc_EMIT(): assert coordinate_vars[1].GetName() == "lat" resampled_ar = ar.GetResampled( - [None, None, ar.GetDimensions()[2]], gdal.GRIORA_NearestNeighbour, None + [None, None, ar.GetDimensions()[2]], + gdal.GRIORA_NearestNeighbour, + None, + ["EMIT_ORTHORECTIFICATION=NO"], ) assert resampled_ar is not None dims = resampled_ar.GetDimensions() @@ -3602,7 +3605,10 @@ def test_netcdf_multidim_getresampled_with_geoloc_EMIT(): assert dims[2].GetSize() == 2 resampled_ar = ar.GetResampled( - [None] * ar.GetDimensionCount(), gdal.GRIORA_NearestNeighbour, None + [None] * ar.GetDimensionCount(), + gdal.GRIORA_NearestNeighbour, + None, + ["EMIT_ORTHORECTIFICATION=NO"], ) assert resampled_ar is not None dims = resampled_ar.GetDimensions() @@ -3646,3 +3652,91 @@ def test_netcdf_multidim_getresampled_with_geoloc_EMIT(): ) == resampled_ar_transposed.Read( array_start_idx=[1, yoff, xoff], count=[1, ysize, xsize] ) + + # Use glt_x and glt_y arrays + resampled_ar = ar.GetResampled( + [None, None, None], gdal.GRIORA_NearestNeighbour, None + ) + assert resampled_ar is not None + dims = resampled_ar.GetDimensions() + assert dims[0].GetName() == "lat" + assert dims[0].GetSize() == 3 + assert dims[1].GetName() == "lon" + assert dims[1].GetSize() == 3 + assert dims[2].GetName() == "bands" + assert dims[2].GetSize() == 2 + assert resampled_ar.GetDataType() == ar.GetDataType() + assert resampled_ar.GetBlockSize() == [3, 3, 2] + assert resampled_ar.GetSpatialRef().GetAuthorityCode(None) == "4326" + assert resampled_ar.GetNoDataValue() == ar.GetNoDataValue() + assert resampled_ar.GetUnit() == ar.GetUnit() + assert ( + resampled_ar.GetAttribute("long_name").ReadAsString() + == ar.GetAttribute("long_name").ReadAsString() + ) + assert struct.unpack("f" * (3 * 3 * 2), resampled_ar.Read()) == ( + -9999.0, + -9999.0, + -9999.0, + -9999.0, + -9999.0, + -9999.0, + -9999.0, + -9999.0, + 30.0, + -30.0, + 40.0, + -40.0, + -9999.0, + -9999.0, + 10.0, + -10.0, + 20.0, + -20.0, + ) + assert struct.unpack( + "d" * (3 * 3 * 2), + resampled_ar.Read( + buffer_datatype=gdal.ExtendedDataType.Create(gdal.GDT_Float64) + ), + ) == ( + -9999.0, + -9999.0, + -9999.0, + -9999.0, + -9999.0, + -9999.0, + -9999.0, + -9999.0, + 30.0, + -30.0, + 40.0, + -40.0, + -9999.0, + -9999.0, + 10.0, + -10.0, + 20.0, + -20.0, + ) + + resampled_ds = resampled_ar.AsClassicDataset(1, 0) + assert resampled_ds.GetGeoTransform() == pytest.approx( + rg.GetAttribute("geotransform").ReadAsDoubleArray() + ) + + assert ( + struct.unpack( + "f", resampled_ar.Read(array_start_idx=[0, 0, 0], count=[1, 1, 1]) + )[0] + == -9999 + ) + assert struct.unpack( + "f" * 4, resampled_ar.Read(array_start_idx=[0, 0, 0], count=[1, 2, 2]) + ) == (-9999, -9999, -9999, -9999) + assert ( + struct.unpack( + "f", resampled_ar.Read(array_start_idx=[2, 1, 1], count=[1, 1, 1]) + )[0] + == -10 + ) diff --git a/frmts/netcdf/netcdfmultidim.cpp b/frmts/netcdf/netcdfmultidim.cpp index f2c0e10c2334..22a5316d5487 100644 --- a/frmts/netcdf/netcdfmultidim.cpp +++ b/frmts/netcdf/netcdfmultidim.cpp @@ -684,6 +684,12 @@ class netCDFVariable final : public GDALPamMDArray, public netCDFAttributeHolder } bool Rename(const std::string &osNewName) override; + + std::shared_ptr + GetResampled(const std::vector> &apoNewDims, + GDALRIOResampleAlg resampleAlg, + const OGRSpatialReference *poTargetSRS, + CSLConstList papszOptions) const override; }; /************************************************************************/ @@ -4211,6 +4217,427 @@ void netCDFVariable::NotifyChildrenOfRenaming() iter.second->ParentRenamed(m_osFullName); } +/************************************************************************/ +/* GetPAM() */ +/************************************************************************/ + +static std::shared_ptr +GetPAM(const std::shared_ptr &poParent) +{ + auto poPamArray = dynamic_cast(poParent.get()); + if (poPamArray) + return poPamArray->GetPAM(); + return nullptr; +} + +/************************************************************************/ +/* GetResampled() */ +/************************************************************************/ + +class GLTOrthoRectifiedArray final : public GDALPamMDArray +{ + private: + std::shared_ptr m_poParent{}; + std::vector> m_apoDims; + std::vector m_anBlockSize; + GDALExtendedDataType m_dt; + std::shared_ptr m_poSRS{}; + std::shared_ptr m_poVarX{}; + std::shared_ptr m_poVarY{}; + std::shared_ptr m_poGLTX{}; + std::shared_ptr m_poGLTY{}; + int m_nGLTIndexOffset = 0; + + protected: + GLTOrthoRectifiedArray( + const std::shared_ptr &poParent, + const std::vector> &apoDims, + const std::vector &anBlockSize) + : GDALAbstractMDArray(std::string(), "GLTOrthoRectifiedArray view of " + + poParent->GetFullName()), + GDALPamMDArray(std::string(), + "GLTOrthoRectifiedArray view of " + + poParent->GetFullName(), + ::GetPAM(poParent)), + m_poParent(std::move(poParent)), m_apoDims(apoDims), + m_anBlockSize(anBlockSize), m_dt(m_poParent->GetDataType()) + { + CPLAssert(apoDims.size() == m_poParent->GetDimensionCount()); + CPLAssert(anBlockSize.size() == m_poParent->GetDimensionCount()); + } + + bool IRead(const GUInt64 *arrayStartIdx, const size_t *count, + const GInt64 *arrayStep, const GPtrDiff_t *bufferStride, + const GDALExtendedDataType &bufferDataType, + void *pDstBuffer) const override; + + public: + static std::shared_ptr + Create(const std::shared_ptr &poParent, + const std::shared_ptr &poGLTX, + const std::shared_ptr &poGLTY, int nGLTIndexOffset, + const std::vector &adfGeoTransform) + { + std::vector> apoNewDims; + + auto poDimY = std::make_shared( + std::string(), CF_LATITUDE_VAR_NAME, GDAL_DIM_TYPE_HORIZONTAL_Y, + "NORTH", poGLTX->GetDimensions()[0]->GetSize()); + auto varY = GDALMDArrayRegularlySpaced::Create( + std::string(), poDimY->GetName(), poDimY, + adfGeoTransform[3] + adfGeoTransform[5] / 2, adfGeoTransform[5], 0); + poDimY->SetIndexingVariable(varY); + apoNewDims.emplace_back(poDimY); + + auto poDimX = std::make_shared( + std::string(), CF_LONGITUDE_VAR_NAME, GDAL_DIM_TYPE_HORIZONTAL_X, + "EAST", poGLTX->GetDimensions()[1]->GetSize()); + auto varX = GDALMDArrayRegularlySpaced::Create( + std::string(), poDimX->GetName(), poDimX, + adfGeoTransform[0] + adfGeoTransform[1] / 2, adfGeoTransform[1], 0); + poDimX->SetIndexingVariable(varX); + apoNewDims.emplace_back(poDimX); + + apoNewDims.emplace_back(poParent->GetDimensions()[2]); + + const std::vector anBlockSize = std::vector{ + std::min(apoNewDims[0]->GetSize(), 512), + std::min(apoNewDims[1]->GetSize(), 512), + poParent->GetDimensions()[2]->GetSize()}; + + auto newAr(std::shared_ptr( + new GLTOrthoRectifiedArray(poParent, apoNewDims, anBlockSize))); + newAr->SetSelf(newAr); + newAr->m_poVarX = varX; + newAr->m_poVarY = varY; + newAr->m_poGLTX = poGLTX; + newAr->m_poGLTY = poGLTY; + newAr->m_nGLTIndexOffset = nGLTIndexOffset; + OGRSpatialReference oSRS; + oSRS.importFromEPSG(4326); + newAr->m_poSRS.reset(oSRS.Clone()); + newAr->m_poSRS->SetDataAxisToSRSAxisMapping( + {/*latIdx = */ 0, /* lonIdx = */ 1}); + return newAr; + } + + bool IsWritable() const override + { + return false; + } + + const std::string &GetFilename() const override + { + return m_poParent->GetFilename(); + } + + const std::vector> & + GetDimensions() const override + { + return m_apoDims; + } + + const GDALExtendedDataType &GetDataType() const override + { + return m_dt; + } + + std::shared_ptr GetSpatialRef() const override + { + return m_poSRS; + } + + std::vector GetBlockSize() const override + { + return m_anBlockSize; + } + + std::shared_ptr + GetAttribute(const std::string &osName) const override + { + return m_poParent->GetAttribute(osName); + } + + std::vector> + GetAttributes(CSLConstList papszOptions = nullptr) const override + { + return m_poParent->GetAttributes(papszOptions); + } + + const std::string &GetUnit() const override + { + return m_poParent->GetUnit(); + } + + const void *GetRawNoDataValue() const override + { + return m_poParent->GetRawNoDataValue(); + } + + double GetOffset(bool *pbHasOffset, + GDALDataType *peStorageType) const override + { + return m_poParent->GetOffset(pbHasOffset, peStorageType); + } + + double GetScale(bool *pbHasScale, + GDALDataType *peStorageType) const override + { + return m_poParent->GetScale(pbHasScale, peStorageType); + } +}; + +/************************************************************************/ +/* GDALMDArrayResampled::IRead() */ +/************************************************************************/ + +bool GLTOrthoRectifiedArray::IRead(const GUInt64 *arrayStartIdx, + const size_t *count, const GInt64 *arrayStep, + const GPtrDiff_t *bufferStride, + const GDALExtendedDataType &bufferDataType, + void *pDstBuffer) const +{ + if (bufferDataType.GetClass() != GEDTC_NUMERIC) + return false; + + const size_t nXYValsCount = count[0] * count[1]; + const auto eInt32DT = GDALExtendedDataType::Create(GDT_Int32); + std::vector anGLTX; + std::vector anGLTY; + try + { + anGLTX.resize(nXYValsCount); + anGLTY.resize(nXYValsCount); + } + catch (const std::bad_alloc &e) + { + CPLError(CE_Failure, CPLE_OutOfMemory, + "GLTOrthoRectifiedArray::IRead(): %s", e.what()); + return false; + } + if (!m_poGLTX->Read(arrayStartIdx, count, arrayStep, nullptr, eInt32DT, + anGLTX.data()) || + !m_poGLTY->Read(arrayStartIdx, count, arrayStep, nullptr, eInt32DT, + anGLTY.data())) + { + return false; + } + + int32_t nMinX = std::numeric_limits::max(); + int32_t nMaxX = std::numeric_limits::min(); + const auto nXSize = m_poParent->GetDimensions()[0]->GetSize(); + for (size_t i = 0; i < nXYValsCount; ++i) + { + const int64_t nX64 = + static_cast(anGLTX[i]) + m_nGLTIndexOffset; + if (nX64 >= 0 && static_cast(nX64) < nXSize) + { + const int32_t nX = static_cast(nX64); + if (nX < nMinX) + nMinX = nX; + if (nX > nMaxX) + nMaxX = nX; + } + } + + int32_t nMinY = std::numeric_limits::max(); + int32_t nMaxY = std::numeric_limits::min(); + const auto nYSize = m_poParent->GetDimensions()[0]->GetSize(); + for (size_t i = 0; i < nXYValsCount; ++i) + { + const int64_t nY64 = + static_cast(anGLTY[i]) + m_nGLTIndexOffset; + if (nY64 >= 0 && static_cast(nY64) < nYSize) + { + const int32_t nY = static_cast(nY64); + if (nY < nMinY) + nMinY = nY; + if (nY > nMaxY) + nMaxY = nY; + } + } + + const auto eBufferDT = bufferDataType.GetNumericDataType(); + auto pRawNoDataValue = GetRawNoDataValue(); + std::vector abyNoData(16); + if (pRawNoDataValue) + GDALCopyWords(pRawNoDataValue, GetDataType().GetNumericDataType(), 0, + abyNoData.data(), eBufferDT, 0, 1); + + const auto nBufferDTSize = bufferDataType.GetSize(); + if (nMinX > nMaxX || nMinY > nMaxY) + { + for (size_t iY = 0; iY < count[0]; ++iY) + { + for (size_t iX = 0; iX < count[1]; ++iX) + { + GByte *pabyDstBuffer = + static_cast(pDstBuffer) + + (iY * bufferStride[0] + iX * bufferStride[1]) * + static_cast(nBufferDTSize); + GDALCopyWords(abyNoData.data(), eBufferDT, 0, pabyDstBuffer, + eBufferDT, + static_cast(bufferStride[2] * nBufferDTSize), + static_cast(count[2])); + } + } + return true; + } + + GUInt64 parentArrayIdxStart[3] = {static_cast(nMinY), + static_cast(nMinX), + arrayStartIdx[2]}; + size_t parentCount[3] = {static_cast(nMaxY - nMinY + 1), + static_cast(nMaxX - nMinX + 1), count[2]}; + GInt64 parentArrayStep[3] = {1, 1, arrayStep[2]}; + + size_t nParentValueSize = nBufferDTSize; + for (int i = 0; i < 3; ++i) + { + if (parentCount[i] > + std::numeric_limits::max() / nParentValueSize) + { + CPLError( + CE_Failure, CPLE_OutOfMemory, + "GLTOrthoRectifiedArray::IRead(): too big temporary array"); + return false; + } + nParentValueSize *= parentCount[i]; + } + + GPtrDiff_t parentStride[3] = { + static_cast(parentCount[1] * parentCount[2]), + static_cast(parentCount[2]), 1}; + std::vector parentValues; + try + { + parentValues.resize(nParentValueSize); + } + catch (const std::bad_alloc &e) + { + CPLError(CE_Failure, CPLE_OutOfMemory, + "GLTOrthoRectifiedArray::IRead(): %s", e.what()); + return false; + } + if (!m_poParent->Read(parentArrayIdxStart, parentCount, parentArrayStep, + parentStride, bufferDataType, parentValues.data())) + { + return false; + } + + size_t iGLTIndex = 0; + for (size_t iY = 0; iY < count[0]; ++iY) + { + for (size_t iX = 0; iX < count[1]; ++iX, ++iGLTIndex) + { + const int64_t nX64 = + static_cast(anGLTX[iGLTIndex]) + m_nGLTIndexOffset; + const int64_t nY64 = + static_cast(anGLTY[iGLTIndex]) + m_nGLTIndexOffset; + GByte *pabyDstBuffer = + static_cast(pDstBuffer) + + (iY * bufferStride[0] + iX * bufferStride[1]) * + static_cast(nBufferDTSize); + if (nX64 >= nMinX && nX64 <= nMaxX && nY64 >= nMinY && + nY64 <= nMaxY) + { + const int32_t iSrcX = static_cast(nX64) - nMinX; + const int32_t iSrcY = static_cast(nY64) - nMinY; + const GByte *pabySrcBuffer = + parentValues.data() + (iSrcY * parentCount[1] + iSrcX) * + parentCount[2] * nBufferDTSize; + GDALCopyWords(pabySrcBuffer, eBufferDT, + static_cast(nBufferDTSize), pabyDstBuffer, + eBufferDT, + static_cast(bufferStride[2] * nBufferDTSize), + static_cast(count[2])); + } + else + { + GDALCopyWords(abyNoData.data(), eBufferDT, 0, pabyDstBuffer, + eBufferDT, + static_cast(bufferStride[2] * nBufferDTSize), + static_cast(count[2])); + } + } + } + + return true; +} + +/************************************************************************/ +/* GetResampled() */ +/************************************************************************/ + +std::shared_ptr netCDFVariable::GetResampled( + const std::vector> &apoNewDims, + GDALRIOResampleAlg resampleAlg, const OGRSpatialReference *poTargetSRS, + CSLConstList papszOptions) const +{ + // Special case for NASA EMIT datasets + auto apoDims = GetDimensions(); + if (apoDims.size() == 3 && apoDims[0]->GetName() == "downtrack" && + apoDims[1]->GetName() == "crosstrack" && + apoDims[2]->GetName() == "bands" && poTargetSRS == nullptr && + (apoNewDims == std::vector>(3) || + apoNewDims == + std::vector>{nullptr, nullptr, + apoDims[2]}) && + CPLTestBool(CSLFetchNameValueDef(papszOptions, + "EMIT_ORTHORECTIFICATION", "YES"))) + { + auto poRootGroup = netCDFGroup::Create(m_poShared, nullptr, m_gid); + if (poRootGroup) + { + auto poAttrGeotransform = poRootGroup->GetAttribute("geotransform"); + auto poLocationGroup = poRootGroup->OpenGroup("location"); + if (poAttrGeotransform && + poAttrGeotransform->GetDataType().GetClass() == GEDTC_NUMERIC && + poAttrGeotransform->GetDimensionCount() == 1 && + poAttrGeotransform->GetDimensionsSize()[0] == 6 && + poLocationGroup) + { + auto poGLT_X = poLocationGroup->OpenMDArray("glt_x"); + auto poGLT_Y = poLocationGroup->OpenMDArray("glt_y"); + if (poGLT_X && poGLT_X->GetDimensionCount() == 2 && + poGLT_X->GetDimensions()[0]->GetName() == "ortho_y" && + poGLT_X->GetDimensions()[1]->GetName() == "ortho_x" && + poGLT_Y && poGLT_Y->GetDimensionCount() == 2 && + poGLT_Y->GetDimensions()[0]->GetName() == "ortho_y" && + poGLT_Y->GetDimensions()[1]->GetName() == "ortho_x") + { + auto self = + std::dynamic_pointer_cast(m_pSelf.lock()); + CPLAssert(self); + if (GetDataType().GetClass() != GEDTC_NUMERIC) + { + CPLError( + CE_Failure, CPLE_AppDefined, + "GetResampled() only supports numeric data type"); + return nullptr; + } + return GLTOrthoRectifiedArray::Create( + self, poGLT_X, poGLT_Y, + /* nGLTIndexOffset = */ -1, + poAttrGeotransform->ReadAsDoubleArray()); + } + } + } + } + + if (CPLTestBool(CSLFetchNameValueDef(papszOptions, + "EMIT_ORTHORECTIFICATION", "NO"))) + { + CPLError(CE_Failure, CPLE_AppDefined, + "EMIT_ORTHORECTIFICATION required, but dataset and/or " + "parameters are not compatible with it"); + return nullptr; + } + + return GDALMDArray::GetResampled(apoNewDims, resampleAlg, poTargetSRS, + papszOptions); +} + /************************************************************************/ /* retrieveAttributeParentName() */ /************************************************************************/ diff --git a/gcore/gdal_priv.h b/gcore/gdal_priv.h index 01e894b2dba7..7a21c9415ba3 100644 --- a/gcore/gdal_priv.h +++ b/gcore/gdal_priv.h @@ -2973,7 +2973,7 @@ class CPL_DLL GDALMDArray : virtual public GDALAbstractMDArray, virtual std::shared_ptr GetMask(CSLConstList papszOptions) const; - std::shared_ptr + virtual std::shared_ptr GetResampled(const std::vector> &apoNewDims, GDALRIOResampleAlg resampleAlg, const OGRSpatialReference *poTargetSRS, diff --git a/gcore/gdalmultidim.cpp b/gcore/gdalmultidim.cpp index ec37cd92458e..fd727c808b75 100644 --- a/gcore/gdalmultidim.cpp +++ b/gcore/gdalmultidim.cpp @@ -8195,6 +8195,9 @@ bool GDALMDArrayResampled::IRead(const GUInt64 *arrayStartIdx, * Currently this method can only resample along the last 2 dimensions, unless * orthorectifying a NASA EMIT dataset. * + * For NASA EMIT netCDF datasets, if apoNewDims[] and poTargetSRS is NULL, the + * geometry lookup table (GLT) is used for fast orthorectification. + * * @param apoNewDims New dimensions. Its size should be GetDimensionCount(). * apoNewDims[i] can be NULL to let the method automatically * determine it. From 337aab7f3b73a746874e03291f80e62fa546c968 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 20 Sep 2023 12:39:57 +0200 Subject: [PATCH 5/6] gdalmdimtranslate: add support for resample=yes array spec option --- apps/gdalmdimtranslate_lib.cpp | 105 +++++++++++++----- .../utilities/test_gdalmdimtranslate_lib.py | 49 ++++++++ doc/source/programs/gdalmdimtranslate.rst | 12 +- frmts/vrt/vrtmultidim.cpp | 23 +++- 4 files changed, 154 insertions(+), 35 deletions(-) diff --git a/apps/gdalmdimtranslate_lib.cpp b/apps/gdalmdimtranslate_lib.cpp index e8e87d093fd8..ec22d513c3f2 100644 --- a/apps/gdalmdimtranslate_lib.cpp +++ b/apps/gdalmdimtranslate_lib.cpp @@ -602,7 +602,7 @@ static bool ParseArraySpec(const std::string &arraySpec, std::string &srcName, std::string &dstName, int &band, std::vector &anTransposedAxis, std::string &viewExpr, - GDALExtendedDataType &outputType) + GDALExtendedDataType &outputType, bool &bResampled) { if (!STARTS_WITH(arraySpec.c_str(), "name=") && !STARTS_WITH(arraySpec.c_str(), "band=")) @@ -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, @@ -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 &poSrcArrayIn, + const std::string &arraySpec, const std::shared_ptr &poSrcRootGroup, const std::shared_ptr &poSrcGroup, const std::shared_ptr &poDstRootGroup, @@ -724,14 +730,16 @@ static bool TranslateArray( int band = -1; std::vector 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 srcArray; + bool bSrcArrayAccessibleThroughSrcGroup = true; if (poSrcRootGroup && poSrcGroup) { if (!srcArrayName.empty() && srcArrayName[0] == '/') @@ -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 @@ -754,6 +770,17 @@ static bool TranslateArray( } auto tmpArray = srcArray; + + if (bResampled) + { + tmpArray = + tmpArray->GetResampled(std::vector>( + tmpArray->GetDimensionCount()), + GRIORA_NearestNeighbour, nullptr, nullptr); + if (!tmpArray) + return false; + } + if (!anTransposedAxis.empty()) { tmpArray = tmpArray->Transpose(anTransposedAxis); @@ -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)) @@ -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) { @@ -1139,23 +1168,39 @@ static bool TranslateArray( } } - const auto dimCount(tmpArray->GetDimensionCount()); - std::vector anSrcOffset(dimCount); - std::vector 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( + dfStart, dfIncrement); + dstArrayVRT->AddSource(std::move(poSource)); + } + else { - anCount[i] = tmpArrayDims[i]->GetSize(); + const auto dimCount(tmpArray->GetDimensionCount()); + std::vector anSrcOffset(dimCount); + std::vector anCount(dimCount); + for (size_t i = 0; i < dimCount; ++i) + { + anCount[i] = tmpArrayDims[i]->GetSize(); + } + std::vector anStep(dimCount, 1); + std::vector anDstOffset(dimCount); + std::unique_ptr 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 anStep(dimCount, 1); - std::vector anDstOffset(dimCount); - std::unique_ptr 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; } @@ -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; @@ -1431,10 +1476,10 @@ static bool TranslateInternal(std::shared_ptr &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; } @@ -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 diff --git a/autotest/utilities/test_gdalmdimtranslate_lib.py b/autotest/utilities/test_gdalmdimtranslate_lib.py index 8f5e9591616d..e7b2e00e17d0 100755 --- a/autotest/utilities/test_gdalmdimtranslate_lib.py +++ b/autotest/utilities/test_gdalmdimtranslate_lib.py @@ -871,6 +871,55 @@ def test_gdalmdimtranslate_array_with_view(): assert dims[1].GetSize() == 5 +@pytest.mark.require_driver("netCDF") +def test_gdalmdimtranslate_array_resample(): + ds = gdal.MultiDimTranslate( + "", + "../gdrivers/data/netcdf/fake_EMIT.nc", + arraySpecs=["name=reflectance,resample=true"], + format="MEM", + ) + rg = ds.GetRootGroup() + resampled_ar = rg.OpenMDArray("reflectance") + dims = resampled_ar.GetDimensions() + assert dims[0].GetName() == "lat" + assert dims[0].GetSize() == 3 + assert dims[1].GetName() == "lon" + assert dims[1].GetSize() == 3 + assert dims[2].GetName() == "bands" + assert dims[2].GetSize() == 2 + assert resampled_ar.GetDataType() == gdal.ExtendedDataType.Create(gdal.GDT_Float32) + assert resampled_ar.GetSpatialRef().GetAuthorityCode(None) == "4326" + assert struct.unpack("f" * (3 * 3 * 2), resampled_ar.Read()) == ( + -9999.0, + -9999.0, + -9999.0, + -9999.0, + -9999.0, + -9999.0, + -9999.0, + -9999.0, + 30.0, + -30.0, + 40.0, + -40.0, + -9999.0, + -9999.0, + 10.0, + -10.0, + 20.0, + -20.0, + ) + + lat = dims[0].GetIndexingVariable() + assert lat + assert struct.unpack("d" * 3, lat.Read()) == (3.5, 2.5, 1.5) + + lon = dims[1].GetIndexingVariable() + assert lon + assert struct.unpack("d" * 3, lon.Read()) == (1.5, 2.5, 3.5) + + def XXXX_test_all(): while True: test_gdalmdimtranslate_no_arg() diff --git a/doc/source/programs/gdalmdimtranslate.rst b/doc/source/programs/gdalmdimtranslate.rst index 0107ab3bf10a..a351844bab5e 100644 --- a/doc/source/programs/gdalmdimtranslate.rst +++ b/doc/source/programs/gdalmdimtranslate.rst @@ -77,12 +77,16 @@ The following command line parameters can appear in any order. may be just an array name, potentially using a fully qualified syntax (/group/subgroup/array_name). Or it can be a combination of options with the syntax: - name={src_array_name}[,dstname={dst_array_name}][,transpose=[{axis1},{axis2},...][,view={view_expr}] + name={src_array_name}[,dstname={dst_array_name}][,resample=yes][,transpose=[{axis1},{axis2},...][,view={view_expr}] - [{axis1},{axis2},...] is the argument of :cpp:func:`GDALMDArray::Transpose`. - For example, transpose=[1,0] switches the axis order of a 2D array. + The following options are processed in that order: - {view_expr} is the value of the *viewExpr* argument of :cpp:func:`GDALMDArray::GetView` + - ``resample=yes`` asks for the array to run through :cpp:func:`GDALMDArray::GetResampled`. + + - [{axis1},{axis2},...] is the argument of :cpp:func:`GDALMDArray::Transpose`. + For example, transpose=[1,0] switches the axis order of a 2D array. + + - {view_expr} is the value of the *viewExpr* argument of :cpp:func:`GDALMDArray::GetView` When specifying a view_expr that performs a slicing or subsetting on a dimension, the equivalent operation will be applied to the corresponding indexing variable. diff --git a/frmts/vrt/vrtmultidim.cpp b/frmts/vrt/vrtmultidim.cpp index 50549835a25c..71eec08b0bdc 100644 --- a/frmts/vrt/vrtmultidim.cpp +++ b/frmts/vrt/vrtmultidim.cpp @@ -1976,6 +1976,25 @@ bool VRTMDArraySourceFromArray::Read(const GUInt64 *arrayStartIdx, poArray = poBand->AsMDArray(); CPLAssert(poArray); } + + std::string osViewExpr = m_osViewExpr; + if (STARTS_WITH(osViewExpr.c_str(), "resample=true,") || + osViewExpr == "resample=true") + { + poArray = + poArray->GetResampled(std::vector>( + poArray->GetDimensionCount()), + GRIORA_NearestNeighbour, nullptr, nullptr); + if (poArray == nullptr) + { + return false; + } + if (osViewExpr == "resample=true") + osViewExpr.clear(); + else + osViewExpr = osViewExpr.substr(strlen("resample=true,")); + } + if (!m_anTransposedAxis.empty()) { poArray = poArray->Transpose(m_anTransposedAxis); @@ -1984,9 +2003,9 @@ bool VRTMDArraySourceFromArray::Read(const GUInt64 *arrayStartIdx, return false; } } - if (!m_osViewExpr.empty()) + if (!osViewExpr.empty()) { - poArray = poArray->GetView(m_osViewExpr); + poArray = poArray->GetView(osViewExpr); if (poArray == nullptr) { return false; From a7ec6b4ced0863e9e216e22a982146ab4a94d821 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 20 Sep 2023 12:45:54 +0200 Subject: [PATCH 6/6] netCDFVariable::IReadWrite(): emit error message from NCGetPutVaraFunc() --- frmts/netcdf/netcdfmultidim.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/frmts/netcdf/netcdfmultidim.cpp b/frmts/netcdf/netcdfmultidim.cpp index 22a5316d5487..1e84dcbf6d49 100644 --- a/frmts/netcdf/netcdfmultidim.cpp +++ b/frmts/netcdf/netcdfmultidim.cpp @@ -3211,7 +3211,10 @@ bool netCDFVariable::IReadWrite( int ret = NCGetPutVaraFunc(m_gid, m_varid, startp.data(), count, buffer); if (ret != NC_NOERR) + { + NCDF_ERR(ret); return false; + } if (bIsRead && (!m_bPerfectDataTypeMatch || bufferDataType.GetNumericDataType() != eDT.GetNumericDataType()))