From 006a75040cd0dbc93b043c2a713d422659ae84be Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 10 Jan 2025 17:18:39 +0100 Subject: [PATCH 1/6] GDALCreateGenImgProjTransformer2(): recognize [SRC_SRS_|DST_SRS_][_DATA_AXIS_TO_SRS_AXIS_MAPPING|_AXIS_MAPPING_STRATEGY] --- alg/gdaltransformer.cpp | 44 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 42 insertions(+), 2 deletions(-) diff --git a/alg/gdaltransformer.cpp b/alg/gdaltransformer.cpp index ba2adc3dc73e..8e5e5e2b02a0 100644 --- a/alg/gdaltransformer.cpp +++ b/alg/gdaltransformer.cpp @@ -1951,10 +1951,47 @@ void *GDALCreateGenImgProjTransformer2(GDALDatasetH hSrcDS, GDALDatasetH hDstDS, const char *pszCO = CSLFetchNameValue(papszOptions, "COORDINATE_OPERATION"); + const auto SetAxisMapping = + [papszOptions](OGRSpatialReference &oSRS, const char *pszPrefix) + { + const char *pszMapping = CSLFetchNameValue( + papszOptions, std::string(pszPrefix) + .append("_DATA_AXIS_TO_SRS_AXIS_MAPPING") + .c_str()); + if (pszMapping) + { + CPLStringList aosTokens(CSLTokenizeString2(pszMapping, ",", 0)); + std::vector anMapping; + for (int i = 0; i < aosTokens.size(); ++i) + anMapping.push_back(atoi(aosTokens[i])); + oSRS.SetDataAxisToSRSAxisMapping(anMapping); + } + else + { + const char *pszStrategy = CSLFetchNameValueDef( + papszOptions, + std::string(pszPrefix).append("_AXIS_MAPPING_STRATEGY").c_str(), + "TRADITIONAL_GIS_ORDER"); + if (EQUAL(pszStrategy, "TRADITIONAL_GIS_ORDER")) + oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); + else if (EQUAL(pszStrategy, "AUTHORITY_COMPLIANT")) + oSRS.SetAxisMappingStrategy(OAMS_AUTHORITY_COMPLIANT); + else + { + CPLError(CE_Warning, CPLE_AppDefined, + "Unrecognized value '%s' for %s", pszStrategy, + std::string(pszPrefix) + .append("_AXIS_MAPPING_STRATEGY") + .c_str()); + return false; + } + } + return true; + }; + OGRSpatialReference oSrcSRS; if (pszSrcSRS) { - oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); if (pszSrcSRS[0] != '\0' && oSrcSRS.SetFromUserInput(pszSrcSRS) != OGRERR_NONE) { @@ -1962,12 +1999,13 @@ void *GDALCreateGenImgProjTransformer2(GDALDatasetH hSrcDS, GDALDatasetH hDstDS, "Failed to import coordinate system `%s'.", pszSrcSRS); return nullptr; } + if (!SetAxisMapping(oSrcSRS, "SRC_SRS")) + return nullptr; } OGRSpatialReference oDstSRS; if (pszDstSRS) { - oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); if (pszDstSRS[0] != '\0' && oDstSRS.SetFromUserInput(pszDstSRS) != OGRERR_NONE) { @@ -1975,6 +2013,8 @@ void *GDALCreateGenImgProjTransformer2(GDALDatasetH hSrcDS, GDALDatasetH hDstDS, "Failed to import coordinate system `%s'.", pszDstSRS); return nullptr; } + if (!SetAxisMapping(oDstSRS, "DST_SRS")) + return nullptr; } const char *pszSrcGeolocArray = From 419d118c134e704db1440a53b3c1f190e672d76b Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 9 Jan 2025 20:31:59 +0100 Subject: [PATCH 2/6] Add convenience GDALDataset::GeolocationToPixelLine() and GDALRasterBand::InterpolateAtGeolocation() and corresponding C functions GDALDatasetGeolocationToPixelLine() and GDALRasterInterpolateAtGeolocation() and Python Band.InterpolateAtGeolocation() --- autotest/gcore/interpolateatpoint.py | 82 ++++++++++++++++- doc/source/spelling_wordlist.txt | 6 ++ gcore/gdal.h | 10 +++ gcore/gdal_priv.h | 11 +++ gcore/gdaldataset.cpp | 127 +++++++++++++++++++++++++++ gcore/gdalrasterband.cpp | 94 +++++++++++++++++++- swig/include/Band.i | 23 +++++ swig/include/java/gdal_java.i | 1 + swig/include/python/gdal_python.i | 69 +++++++++++++++ 9 files changed, 420 insertions(+), 3 deletions(-) diff --git a/autotest/gcore/interpolateatpoint.py b/autotest/gcore/interpolateatpoint.py index d3ee6add3249..1409ff6a735e 100755 --- a/autotest/gcore/interpolateatpoint.py +++ b/autotest/gcore/interpolateatpoint.py @@ -15,7 +15,7 @@ import gdaltest import pytest -from osgeo import gdal +from osgeo import gdal, osr ############################################################################### # Test some cases of InterpolateAtPoint @@ -353,3 +353,83 @@ def test_interpolateatpoint_big_complex(): assert res == pytest.approx((182 + 122j), 1e-6) res = ds.GetRasterBand(1).InterpolateAtPoint(255.5, 63.5, gdal.GRIORA_Cubic) assert res == pytest.approx((182 + 122j), 1e-6) + + +############################################################################### +# Test InterpolateAtGeolocation + + +def test_interpolateatgeolocation_1(): + + ds = gdal.Open("data/byte.tif") + + gt = ds.GetGeoTransform() + X = gt[0] + 10 * gt[1] + Y = gt[3] + 12 * gt[5] + + got_nearest = ds.GetRasterBand(1).InterpolateAtGeolocation( + X, Y, None, gdal.GRIORA_NearestNeighbour + ) + assert got_nearest == pytest.approx(173, 1e-6) + got_bilinear = ds.GetRasterBand(1).InterpolateAtGeolocation( + X, Y, None, gdal.GRIORA_Bilinear + ) + assert got_bilinear == pytest.approx(139.75, 1e-6) + got_cubic = ds.GetRasterBand(1).InterpolateAtGeolocation( + X, Y, None, gdal.GRIORA_CubicSpline + ) + assert got_cubic == pytest.approx(138.02, 1e-2) + got_cubic = ds.GetRasterBand(1).InterpolateAtGeolocation( + X, Y, None, gdal.GRIORA_Cubic + ) + assert got_cubic == pytest.approx(145.57, 1e-2) + + wgs84_srs_auth_compliant = osr.SpatialReference() + wgs84_srs_auth_compliant.SetFromUserInput("WGS84") + wgs84_srs_auth_compliant.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT) + ct = osr.CoordinateTransformation(ds.GetSpatialRef(), wgs84_srs_auth_compliant) + wgs84_lat, wgs84_lon, _ = ct.TransformPoint(X, Y, 0) + + got_nearest = ds.GetRasterBand(1).InterpolateAtGeolocation( + wgs84_lat, wgs84_lon, wgs84_srs_auth_compliant, gdal.GRIORA_NearestNeighbour + ) + assert got_nearest == pytest.approx(173, 1e-6) + + wgs84_trad_gis_order = osr.SpatialReference() + wgs84_trad_gis_order.SetFromUserInput("WGS84") + wgs84_trad_gis_order.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) + got_nearest = ds.GetRasterBand(1).InterpolateAtGeolocation( + wgs84_lon, wgs84_lat, wgs84_trad_gis_order, gdal.GRIORA_NearestNeighbour + ) + assert got_nearest == pytest.approx(173, 1e-6) + + wgs84_lon_lat_srs = osr.SpatialReference() + wgs84_lon_lat_srs.SetFromUserInput("WGS84") + wgs84_lon_lat_srs.SetDataAxisToSRSAxisMapping([2, 1]) + got_nearest = ds.GetRasterBand(1).InterpolateAtGeolocation( + wgs84_lon, wgs84_lat, wgs84_lon_lat_srs, gdal.GRIORA_NearestNeighbour + ) + assert got_nearest == pytest.approx(173, 1e-6) + + wgs84_lat_lon_srs = osr.SpatialReference() + wgs84_lat_lon_srs.SetFromUserInput("WGS84") + wgs84_lat_lon_srs.SetDataAxisToSRSAxisMapping([1, 2]) + got_nearest = ds.GetRasterBand(1).InterpolateAtGeolocation( + wgs84_lat, wgs84_lon, wgs84_lat_lon_srs, gdal.GRIORA_NearestNeighbour + ) + assert got_nearest == pytest.approx(173, 1e-6) + + X = gt[0] + -1 * gt[1] + Y = gt[3] + -1 * gt[5] + assert ( + ds.GetRasterBand(1).InterpolateAtGeolocation( + X, Y, None, gdal.GRIORA_NearestNeighbour + ) + is None + ) + + ungeoreferenced_ds = gdal.GetDriverByName("MEM").Create("", 1, 1) + with pytest.raises(Exception, match="Unable to compute a transformation"): + assert ungeoreferenced_ds.GetRasterBand(1).InterpolateAtGeolocation( + 0, 0, None, gdal.GRIORA_NearestNeighbour + ) diff --git a/doc/source/spelling_wordlist.txt b/doc/source/spelling_wordlist.txt index f39a583c6268..b47fdb070d60 100644 --- a/doc/source/spelling_wordlist.txt +++ b/doc/source/spelling_wordlist.txt @@ -639,6 +639,8 @@ dfCutlineBlendDist dfCxtX dfDistance dfErrorThreshold +dfGeolocX +dfGeolocY dfInvisibleVal dfLLX dfMax @@ -1065,6 +1067,8 @@ geoloc geolocated geolocation Geolocation +geolocX +geolocY Géologiques Geomatics Geomedia @@ -2127,6 +2131,7 @@ nYSize nYStride nz nZStride +OAMS oauth OAuth objdir @@ -2910,6 +2915,7 @@ ServerUrl setAccessControl setAccessControlRecursive setargv +SetAxisMappingStrategy setCoordinateDimension setfeature SetFeature diff --git a/gcore/gdal.h b/gcore/gdal.h index 8f83b5d3bc9b..fea2e37f7d10 100644 --- a/gcore/gdal.h +++ b/gcore/gdal.h @@ -1247,6 +1247,10 @@ CPLErr CPL_DLL GDALSetSpatialRef(GDALDatasetH, OGRSpatialReferenceH); CPLErr CPL_DLL CPL_STDCALL GDALGetGeoTransform(GDALDatasetH, double *); CPLErr CPL_DLL CPL_STDCALL GDALSetGeoTransform(GDALDatasetH, double *); +CPLErr CPL_DLL GDALDatasetGeolocationToPixelLine( + GDALDatasetH, double dfGeolocX, double dfGeolocY, OGRSpatialReferenceH hSRS, + double *pdfPixel, double *pdfLine, CSLConstList papszTransformerOptions); + int CPL_DLL CPL_STDCALL GDALGetGCPCount(GDALDatasetH); const char CPL_DLL *CPL_STDCALL GDALGetGCPProjection(GDALDatasetH); OGRSpatialReferenceH CPL_DLL GDALGetGCPSpatialRef(GDALDatasetH); @@ -1703,6 +1707,12 @@ CPLErr CPL_DLL GDALRasterInterpolateAtPoint(GDALRasterBandH hBand, double *pdfRealValue, double *pdfImagValue); +CPLErr CPL_DLL GDALRasterInterpolateAtGeolocation( + GDALRasterBandH hBand, double dfGeolocX, double dfGeolocY, + OGRSpatialReferenceH hSRS, GDALRIOResampleAlg eInterpolation, + double *pdfRealValue, double *pdfImagValue, + CSLConstList papszTransformerOptions); + /** Generic pointer for the working structure of VRTProcessedDataset * function. */ typedef void *VRTPDWorkingDataPtr; diff --git a/gcore/gdal_priv.h b/gcore/gdal_priv.h index 5ec00f099baf..6405ebe9f6e3 100644 --- a/gcore/gdal_priv.h +++ b/gcore/gdal_priv.h @@ -715,6 +715,11 @@ class CPL_DLL GDALDataset : public GDALMajorObject virtual CPLErr GetGeoTransform(double *padfTransform); virtual CPLErr SetGeoTransform(double *padfTransform); + CPLErr GeolocationToPixelLine( + double dfGeolocX, double dfGeolocY, const OGRSpatialReference *poSRS, + double *pdfPixel, double *pdfLine, + CSLConstList papszTransformerOptions = nullptr) const; + virtual CPLErr AddBand(GDALDataType eType, char **papszOptions = nullptr); virtual void *GetInternalHandle(const char *pszHandleName); @@ -1874,6 +1879,12 @@ class CPL_DLL GDALRasterBand : public GDALMajorObject std::shared_ptr AsMDArray() const; + CPLErr InterpolateAtGeolocation( + double dfGeolocX, double dfGeolocY, const OGRSpatialReference *poSRS, + GDALRIOResampleAlg eInterpolation, double *pdfRealValue, + double *pdfImagValue = nullptr, + CSLConstList papszTransfomerOptions = nullptr) const; + virtual CPLErr InterpolateAtPoint(double dfPixel, double dfLine, GDALRIOResampleAlg eInterpolation, double *pdfRealValue, diff --git a/gcore/gdaldataset.cpp b/gcore/gdaldataset.cpp index a6bf53191298..bbb2343e6f0d 100644 --- a/gcore/gdaldataset.cpp +++ b/gcore/gdaldataset.cpp @@ -36,6 +36,7 @@ #include "cpl_string.h" #include "cpl_vsi.h" #include "cpl_vsi_error.h" +#include "gdal_alg.h" #include "ogr_api.h" #include "ogr_attrind.h" #include "ogr_core.h" @@ -10301,3 +10302,129 @@ GDALDataset::Clone(int nScopeFlags, [[maybe_unused]] bool bCanShareState) const } //! @endcond + +/************************************************************************/ +/* GeolocationToPixelLine() */ +/************************************************************************/ + +/** Transform georeferenced coordinates to pixel/line coordinates. + * + * When poSRS is null, those georeferenced coordinates (dfGeolocX, dfGeolocY) + * must be in the "natural" SRS of the dataset, that is the one returned by + * GetSpatialRef() if there is a geotransform, GetGCPSpatialRef() if there are + * GCPs, WGS 84 if there are RPC coefficients, or the SRS of the geolocation + * array (generally WGS 84) if there is a geolocation array. + * If that natural SRS is a geographic one, dfGeolocX must be a longitude, and + * dfGeolocY a latitude. If that natural SRS is a projected one, dfGeolocX must + * be a easting, and dfGeolocY a northing. + * + * When poSRS is set to a non-null value, (dfGeolocX, dfGeolocY) must be + * expressed in that CRS, and that tuple must be conformant with the + * data-axis-to-crs-axis setting of poSRS, that is the one returned by + * the OGRSpatialReference::GetDataAxisToSRSAxisMapping(). If you want to be sure + * of the axis order, then make sure to call poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER) + * before calling this method, and in that case, dfGeolocX must be a longitude + * or an easting value, and dfGeolocX a latitude or a northing value. + * + * This method uses GDALCreateGenImgProjTransformer2() underneath. + * + * @param dfGeolocX X coordinate of the position (longitude or easting if poSRS + * is null, otherwise consistent with poSRS data-axis-to-crs-axis mapping), + * where interpolation should be done. + * @param dfGeolocY Y coordinate of the position (latitude or northing if poSRS + * is null, otherwise consistent with poSRS data-axis-to-crs-axis mapping), + * where interpolation should be done. + * @param poSRS If set, override the natural CRS in which dfGeolocX, dfGeolocY are expressed + * @param[out] pdfPixel Pointer to the variable where to the store the pixel/column coordinate. + * @param[out] pdfLine Pointer to the variable where to the store the line coordinate. + * @param papszTransformerOptions Options accepted by GDALCreateGenImgProjTransformer2(), or nullptr. + * + * @return CE_None on success, or an error code on failure. + * @since GDAL 3.11 + */ + +CPLErr +GDALDataset::GeolocationToPixelLine(double dfGeolocX, double dfGeolocY, + const OGRSpatialReference *poSRS, + double *pdfPixel, double *pdfLine, + CSLConstList papszTransformerOptions) const +{ + CPLStringList aosTO(papszTransformerOptions); + + if (poSRS) + { + const char *const apszOptions[] = {"FORMAT=WKT2", nullptr}; + const std::string osWKT = poSRS->exportToWkt(apszOptions); + aosTO.SetNameValue("DST_SRS", osWKT.c_str()); + const auto eAxisMappingStrategy = poSRS->GetAxisMappingStrategy(); + if (eAxisMappingStrategy == OAMS_TRADITIONAL_GIS_ORDER) + aosTO.SetNameValue("DST_SRS_AXIS_MAPPING_STRATEGY", + "TRADITIONAL_GIS_ORDER"); + else if (eAxisMappingStrategy == OAMS_AUTHORITY_COMPLIANT) + aosTO.SetNameValue("DST_SRS_AXIS_MAPPING_STRATEGY", + "AUTHORITY_COMPLIANT"); + else + { + const auto &anValues = poSRS->GetDataAxisToSRSAxisMapping(); + std::string osVal; + for (int v : anValues) + { + if (!osVal.empty()) + osVal += ','; + osVal += std::to_string(v); + } + aosTO.SetNameValue("DST_SRS_DATA_AXIS_TO_SRS_AXIS_MAPPING", + osVal.c_str()); + } + } + + auto hTransformer = GDALCreateGenImgProjTransformer2( + GDALDataset::ToHandle(const_cast(this)), nullptr, + aosTO.List()); + if (hTransformer == nullptr) + { + return CE_Failure; + } + + double z = 0; + int bSuccess = 0; + GDALGenImgProjTransform(hTransformer, TRUE, 1, &dfGeolocX, &dfGeolocY, &z, + &bSuccess); + GDALDestroyTransformer(hTransformer); + if (bSuccess) + { + if (pdfPixel) + *pdfPixel = dfGeolocX; + if (pdfLine) + *pdfLine = dfGeolocY; + return CE_None; + } + else + { + return CE_Failure; + } +} + +/************************************************************************/ +/* GDALDatasetGeolocationToPixelLine() */ +/************************************************************************/ + +/** Transform georeferenced coordinates to pixel/line coordinates. + * + * @see GDALDataset::GeolocationToPixelLine() + * @since GDAL 3.11 + */ + +CPLErr GDALDatasetGeolocationToPixelLine(GDALDatasetH hDS, double dfGeolocX, + double dfGeolocY, + OGRSpatialReferenceH hSRS, + double *pdfPixel, double *pdfLine, + CSLConstList papszTransformerOptions) +{ + VALIDATE_POINTER1(hDS, "GDALDatasetGeolocationToPixelLine", CE_Failure); + + GDALDataset *poDS = GDALDataset::FromHandle(hDS); + return poDS->GeolocationToPixelLine( + dfGeolocX, dfGeolocY, OGRSpatialReference::FromHandle(hSRS), pdfPixel, + pdfLine, papszTransformerOptions); +} diff --git a/gcore/gdalrasterband.cpp b/gcore/gdalrasterband.cpp index 187f5dddb872..4185b65cab29 100644 --- a/gcore/gdalrasterband.cpp +++ b/gcore/gdalrasterband.cpp @@ -9662,8 +9662,8 @@ std::shared_ptr GDALRasterBand::AsMDArray() const /************************************************************************/ /** - * \brief Interpolates the value between pixels using - * a resampling algorithm + * \brief Interpolates the value between pixels using a resampling algorithm, + * taking pixel/line coordinates as input. * * @param dfPixel pixel coordinate as a double, where interpolation should be done. * @param dfLine line coordinate as a double, where interpolation should be done. @@ -9726,3 +9726,93 @@ CPLErr GDALRasterInterpolateAtPoint(GDALRasterBandH hBand, double dfPixel, return poBand->InterpolateAtPoint(dfPixel, dfLine, eInterpolation, pdfRealValue, pdfImagValue); } + +/************************************************************************/ +/* InterpolateAtGeolocation() */ +/************************************************************************/ + +/** + * \brief Interpolates the value between pixels using a resampling algorithm, + * taking georeferenced coordinates as input. + * + * When poSRS is null, those georeferenced coordinates (dfGeolocX, dfGeolocY) + * must be in the "natural" SRS of the dataset, that is the one returned by + * GetSpatialRef() if there is a geotransform, GetGCPSpatialRef() if there are + * GCPs, WGS 84 if there are RPC coefficients, or the SRS of the geolocation + * array (generally WGS 84) if there is a geolocation array. + * If that natural SRS is a geographic one, dfGeolocX must be a longitude, and + * dfGeolocY a latitude. If that natural SRS is a projected one, dfGeolocX must + * be a easting, and dfGeolocY a northing. + * + * When poSRS is set to a non-null value, (dfGeolocX, dfGeolocY) must be + * expressed in that CRS, and that tuple must be conformant with the + * data-axis-to-crs-axis setting of poSRS, that is the one returned by + * the OGRSpatialReference::GetDataAxisToSRSAxisMapping(). If you want to be sure + * of the axis order, then make sure to call poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER) + * before calling this method, and in that case, dfGeolocX must be a longitude + * or an easting value, and dfGeolocX a latitude or a northing value. + * + * The GDALDataset::GeolocationToPixelLine() will be used to transform from + * (dfGeolocX,dfGeolocY) georeferenced coordinates to (pixel, line). Refer to + * it for details on how that transformation is done. + * + * @param dfGeolocX X coordinate of the position (longitude or easting if poSRS + * is null, otherwise consistent with poSRS data-axis-to-crs-axis mapping), + * where interpolation should be done. + * @param dfGeolocY Y coordinate of the position (latitude or northing if poSRS + * is null, otherwise consistent with poSRS data-axis-to-crs-axis mapping), + * where interpolation should be done. + * @param poSRS If set, override the natural CRS in which dfGeolocX, dfGeolocY are expressed + * @param eInterpolation interpolation type. Only near, bilinear, cubic and cubicspline are allowed. + * @param pdfRealValue pointer to real part of interpolated value + * @param pdfImagValue pointer to imaginary part of interpolated value (may be null if not needed) + * @param papszTransformerOptions Options accepted by GDALDataset::GeolocationToPixelLine() (GDALCreateGenImgProjTransformer2()), or nullptr. + * + * @return CE_None on success, or an error code on failure. + * @since GDAL 3.11 + */ + +CPLErr GDALRasterBand::InterpolateAtGeolocation( + double dfGeolocX, double dfGeolocY, const OGRSpatialReference *poSRS, + GDALRIOResampleAlg eInterpolation, double *pdfRealValue, + double *pdfImagValue, CSLConstList papszTransformerOptions) const +{ + double dfPixel; + double dfLine; + if (poDS->GeolocationToPixelLine(dfGeolocX, dfGeolocY, poSRS, &dfPixel, + &dfLine, + papszTransformerOptions) != CE_None) + { + return CE_Failure; + } + return InterpolateAtPoint(dfPixel, dfLine, eInterpolation, pdfRealValue, + pdfImagValue); +} + +/************************************************************************/ +/* GDALRasterInterpolateAtGeolocation() */ +/************************************************************************/ + +/** + * \brief Interpolates the value between pixels using a resampling algorithm, + * taking georeferenced coordinates as input. + * + * @see GDALRasterBand::InterpolateAtGeolocation() + * @since GDAL 3.11 + */ + +CPLErr GDALRasterInterpolateAtGeolocation(GDALRasterBandH hBand, + double dfGeolocX, double dfGeolocY, + OGRSpatialReferenceH hSRS, + GDALRIOResampleAlg eInterpolation, + double *pdfRealValue, + double *pdfImagValue, + CSLConstList papszTransformerOptions) +{ + VALIDATE_POINTER1(hBand, "GDALRasterInterpolateAtGeolocation", CE_Failure); + + GDALRasterBand *poBand = GDALRasterBand::FromHandle(hBand); + return poBand->InterpolateAtGeolocation( + dfGeolocX, dfGeolocY, OGRSpatialReference::FromHandle(hSRS), + eInterpolation, pdfRealValue, pdfImagValue, papszTransformerOptions); +} diff --git a/swig/include/Band.i b/swig/include/Band.i index 66d809e166c3..d2eac7f958e7 100644 --- a/swig/include/Band.i +++ b/swig/include/Band.i @@ -668,6 +668,29 @@ CPLErr AdviseRead( int xoff, int yoff, int xsize, int ysize, %clear (CPLErr); #endif + +%apply (double *OUTPUT){double *pdfRealValue, double *pdfImagValue}; +#if !defined(SWIGPYTHON) +%apply (IF_ERROR_RETURN_NONE) { (CPLErr) }; +#endif +%apply (char **options) { char ** transformerOptions }; + CPLErr InterpolateAtGeolocation( double geolocX, double geolocY, + OSRSpatialReferenceShadow* srs, + GDALRIOResampleAlg interpolation, + double *pdfRealValue, + double *pdfImagValue, char** transformerOptions = NULL ) { + if (pdfRealValue) *pdfRealValue = 0; + if (pdfImagValue) *pdfImagValue = 0; + return GDALRasterInterpolateAtGeolocation( self, geolocX, geolocY, + (OGRSpatialReferenceH)srs, interpolation, + pdfRealValue, pdfImagValue, transformerOptions ); + } +#if !defined(SWIGPYTHON) +%clear (CPLErr); +%clear char ** transformerOptions; +#endif + + %apply (double *OUTPUT){double *pdfMin, double *pdfMax}; %apply (int *OUTPUT){int *pnMinX, int *pnMinY}; %apply (int *OUTPUT){int *pnMaxX, int *pnMaxY}; diff --git a/swig/include/java/gdal_java.i b/swig/include/java/gdal_java.i index 9efea1534fb8..70637de4a2ba 100644 --- a/swig/include/java/gdal_java.i +++ b/swig/include/java/gdal_java.i @@ -1122,6 +1122,7 @@ import java.awt.Color; %typemap(javaimports) GDALRasterBandShadow %{ import org.gdal.gdalconst.gdalconstConstants; +import org.gdal.osr.SpatialReference; %} diff --git a/swig/include/python/gdal_python.i b/swig/include/python/gdal_python.i index 06eb319fa67c..bfd4a429dd90 100644 --- a/swig/include/python/gdal_python.i +++ b/swig/include/python/gdal_python.i @@ -5024,6 +5024,75 @@ def InterpolateAtPoint(self, *args, **kwargs): return ret[1] %} +%feature("shadow") InterpolateAtGeolocation %{ +def InterpolateAtGeolocation(self, *args, **kwargs): + """Return the interpolated value at georeferenced coordinates. + See :cpp:func:`GDALRasterBand::InterpolateAtGeolocation`. + + When srs is None, those georeferenced coordinates (geolocX, geolocY) + must be in the "natural" SRS of the dataset, that is the one returned by + GetSpatialRef() if there is a geotransform, GetGCPSpatialRef() if there are + GCPs, WGS 84 if there are RPC coefficients, or the SRS of the geolocation + array (generally WGS 84) if there is a geolocation array. + If that natural SRS is a geographic one, geolocX must be a longitude, and + geolocY a latitude. If that natural SRS is a projected one, geolocX must + be a easting, and geolocY a northing. + + When srs is set to a non-None value, (geolocX, geolocY) must be + expressed in that CRS, and that tuple must be conformant with the + data-axis-to-crs-axis setting of srs, that is the one returned by + the :py:func:`osgeo.osr.SpatialReference.GetDataAxisToSRSAxisMapping(). + If you want to be sure of the axis order, then make sure to call + ``srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)`` + before calling this method, and in that case, geolocX must be a longitude + or an easting value, and geolocX a latitude or a northing value. + + Parameters + ---------- + geolocX : float + X coordinate of the position where interpolation should be done. + Longitude or easting in "natural" CRS if `srs` is None, + otherwise consistent with first axis of `srs`, + taking into account the data-axis-to-crs-axis mapping + geolocY : float + Y coordinate of the position where interpolation should be done. + Latitude or northing in "natural" CRS if `srs` is None, + otherwise consistent with second axis of `srs`, + taking into account the data-axis-to-crs-axis mapping + srs : osgeo.osr.SpatialReference + If set, override the natural CRS in which geolocX, geolocY are expressed + interpolation : GRIOResampleAlg (nearest, bilinear, cubic, cubicspline) + + Returns + ------- + float: + Interpolated value, or ``None`` if it has any error. + + Example + ------- + + >>> from osgeo import gdal, osr + >>> with gdal.Open("my.tif") as ds: + ... wgs84_srs = osr.SpatialReference() + ... wgs84_srs.SetFromUserInput("WGS84") + ... wgs84_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) + ... val = ds.GetRasterBand(1).InterpolateAtGeolocation(longitude_degree, + latitude_degree, + wgs84_srs, + gdal.GRIORA_Bilinear) + """ + + ret = $action(self, *args, **kwargs) + if ret[0] != CE_None: + return None + + from . import gdal + if gdal.DataTypeIsComplex(self.DataType): + return complex(ret[1], ret[2]) + else: + return ret[1] +%} + %feature("shadow") ComputeMinMaxLocation %{ def ComputeMinMaxLocation(self, *args, **kwargs): """Compute the min/max values for a band, and their location. From e734b401bedbc3e4b460df278b0ba2020d59ce5c Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 19 Oct 2024 05:13:15 +0200 Subject: [PATCH 3/6] PDF: remove write support for GEO_ENCODING=OGC_BP for 3.11 Several reasons: - support for SRS has always been partial - and the actual trigger for this pull request: lately the company behind the OGC Best Practice encoding has behaved in a unfriendly way with respect to our community. No need to promote their techs further. --- .../data/pdf/test_pdf_composition.pdf | Bin 8305 -> 7709 bytes .../test_pdf_composition_libpng_1_6_40.pdf | Bin 8281 -> 7685 bytes autotest/gdrivers/pdf.py | 253 +------ doc/source/drivers/raster/pdf.rst | 19 +- frmts/pdf/pdfcreatecopy.cpp | 652 +----------------- frmts/pdf/pdfcreatecopy.h | 5 - frmts/pdf/pdfcreatefromcomposition.cpp | 101 +-- frmts/pdf/pdfcreatefromcomposition.h | 7 - frmts/pdf/pdfdrivercore.cpp | 2 - port/cpl_known_config_options.h | 1 - 10 files changed, 30 insertions(+), 1010 deletions(-) diff --git a/autotest/gdrivers/data/pdf/test_pdf_composition.pdf b/autotest/gdrivers/data/pdf/test_pdf_composition.pdf index a3d82caa93ae4518094f0c3f611b91b5ce2b699e..c59dcdc734b344d4f4fa8a66c3caa4056f84d0e8 100644 GIT binary patch delta 1512 zcmZuxOKTHR6lOvncJeYNYLiB|($qd6ow<*BX?aC5Ds?!>25@}+w8{@hg z7iaBHP%7EFaw90Zb|ENk6-AdK;=MCTnq+bb4Eg4K=X~co=gzPBXER6d`Xv+o4qk#k zgCU5AGd`l$i{WW0PE@3#^45JBAw-BN+oikWAHncth`@I_1-cv!P#%;)%FP{w4e*DO zK0FBn@+=(51KsZ?eT2djc@$3ihhog(5zTpIc>-Po(*6{@4W!^jC;?CXG@!R9|KK0x zN`L!5LXFqvIf6z9xi6T8uO|gL28CcWNUU}?mnLCI1()Dkk9SiN0X?)K34IgrDHw*= zp&l9FUnm5Jp#ThrHzbYQq)6vQ2VR$`kS^ev8TdJLHaxANxw)2wR>Ou9Nh(avAd20G z>t@rAsmQ2gz*=NnGVA-1w6tm-*TRH0(^1i=sI`emhpoglsAmn7g_oydy2zM8*=8y% z7_MxjER;&+Y6US@qOw>nRqc`hFkFnz!(kM2*Iq@ho)+Spq!qigQ@w{6g?LQNYc@9G z!&d!vEYcrk=yF-UB`>Wm%F0T9%ig61JdU4?KqFB+zX+7!bGJdY>)OoUB%A0M73 zj#g!PJ-4`!Bbt`Glh0w^qvJK-!R7TrX}?e@Z0*{JZ)r4PPGAf1pR=Kjf{41NDC-Bs zYsLJIO&OvLDhXw!wN4uc4{W}X{AHzWYD{?BY;)~lv4S=^;FzC5D9a;JQ#oJdABpsc zI}GJI4SL7O zsL$t8HLNjD7`A$0I#6M^v#Lom_%NM{13u$NGrj)l8_EQWu$79*i0~ z-=?lOvxh*7wQ)s#eyD1+fej3M6H;P%icq}=W?()-agLj#rnobm7^@M?FG&~3TlQk! vj+w?hwrZB;v0t+ck7%73P7N4GYQp)#RI2&C>Z3h-ThU{heWb^}w>4#rJP zdSQ%Ldf~a&Utr8eet;%kcwwR!>xE|5Wc>+V=!J0;-!mWdV>%F~ea>^?l)4x*XfyPgM@@I$vv& z!GUryn=Tb{`#p+EHC>|GHZk-hQ8YrRqSEcuVI^BAI-C?V+>B@7IM zfsr{g0>J+cLmMUl#V5cl-I3&PWODn>QyYnIne7e_A16s^i~Tpg%zhge#|UPBMBW`E z($(wJ$g~cOnrMjpW!!nRQU)<=Iz5#?|fFL}mh?ZkFJ8}u2$xvAQP^n6ICh=Tv9Ai!Uv6IXpSxBEZ)TPysh!OQINq8=w0!LMU4 zFCn6?cw+hCMed^BFsXv5!swEizqXloxxr$m&rC9(f#2rA!Zh_&ws-ii?Z1g8q6*(Zks>=rZ<)*x%axZu6dfWc6 RZQd6&iZx+#GhMhX{14{m`#}Hz diff --git a/autotest/gdrivers/data/pdf/test_pdf_composition_libpng_1_6_40.pdf b/autotest/gdrivers/data/pdf/test_pdf_composition_libpng_1_6_40.pdf index 57ae04340d089478391b79aee311f80bc6f7efe1..c4d1b36828026fcd8342bac37790142476b56147 100644 GIT binary patch delta 1479 zcmZuxOH30{6wROne8ZH_lmMl|laB(CblzuvN(hJ*Klr1%B1j63h@nHIqrpVN%9YEz zGI8l%lP2z6kw{`(7>z3>acN9knW%}rH>JRIXq&uw=YGz;=ice(u_q(1wrUjz{??s_ z-*pMl8mDWhQOY;=C>_+m240!JhcL#JSjh8l%Nr%|;RuDVs1Jsa-oS*@G|pUG#KeYw zS{j~d8q^}%LIj9um3gp+GT25ZlRONgW_Zkot`v@`eTps2y|ASvVO@0_tZ4RD_@pKv ztG2)|brXuBgda`y32Ikm(LvQx@V2Qh%rmR$8g-+IWSI%**BtnAfOSZrFqF6fgNX+C zj@rX4uM-DkX=HkRO8&YW|HJ-+dl#jAn zdvfT|-`pf5v6@dPPU%+ju(ILor6rW@Qp^4%+YF`=X?L_akO$9J!g*!XB3u$o@0o}c zS7sg=!I^Qt-CqfcIP@1&pIGq&&zAs!9(@c}bt0sDp`UM+4OOD1Jby8`i#Y=&+1usd z#Ui?oB#nU`OvHDgfHD^f`9+Vor8CLq+98h5uAu8EdwBw7rtDn>POc%Dh^+S}qpkkgfZ?o zc^$1D~Gn48ZhE70)0j&scRdC6F3SZ;3@bnnb z`PY-n{Mp3xoH~lN#($Yu;onClqeS!16Ww1#;~_%%pOZKEH=#NH&Qy}yAvz_x_mYz+ ze;kVN&m(jEY3M1BhHeN`5n>M6nCHLCi~OUo$y4F^X^IDph)5&+-|!TFG8yNO!Z!nq zTN35Z!X=6M$vg@#@K2^j*qNcCbUlHn1~MK=I;glT#YI{|wm}i8!OtQo{tcN-CAF=-Rh7`5XC@L9&dO*|)vhpflpl2X{CKb6 za9lVj_P11_Y;Ss`-MHKET8MSbDVc`rIj#T*|Le+~$ac5QRKI8(&r*=_bT%Iq`Spx! zJ8so`4>87nm|cF+Le-Z1fh@dWWcowtq%JFkdZp?Rd@(jRJAO+;n+^GOIeRxFE8CT0 z=Zu;HfwqL};#BI-(_sA7SnL&g$yxTkm$SW!=lG*6XNw!u$xzAJN161Is#-@U6|o>n z`KPgXwA*kwyXJmaYt@d=90Z;JUCdl17{R~aj=eA-693zNQ1?To8iFoZiYrjHt@ukZ z(c-14PQKAf^k+#nTU(8aSXJ zHE952MAT{UENc@D8e|A$IN*dC493F{&(MQTOfeb;c^b9?JY8dcF#PfL>6k_WrPWQM f@o%rQMO}=mmRD(d_nXeUvO!ELudHNuUz7g_1s(hA diff --git a/autotest/gdrivers/pdf.py b/autotest/gdrivers/pdf.py index 6239ea6b8bf1..2aa26339d5e2 100755 --- a/autotest/gdrivers/pdf.py +++ b/autotest/gdrivers/pdf.py @@ -348,81 +348,6 @@ def test_pdf_iso32000_dpi_300(poppler_or_pdfium): ) -############################################################################### -# Test write support with OGC_BP geo encoding - - -def test_pdf_ogcbp(poppler_or_pdfium_or_podofo): - with gdal.config_option("GDAL_PDF_OGC_BP_WRITE_WKT", "FALSE"): - tst = gdaltest.GDALTest( - "PDF", "byte.tif", 1, None, options=["GEO_ENCODING=OGC_BP"] - ) - gdal.ErrorReset() - tst.testCreateCopy( - check_minmax=0, - check_gt=1, - check_srs=True, - check_checksum_not_null=pdf_checksum_available(), - ) - assert gdal.GetLastErrorMsg() == "" - - -############################################################################### -# Test write support with OGC_BP geo encoding, with DPI=300 - - -def test_pdf_ogcbp_dpi_300(poppler_or_pdfium): - with gdal.config_option("GDAL_PDF_OGC_BP_WRITE_WKT", "FALSE"): - tst = gdaltest.GDALTest( - "PDF", "byte.tif", 1, None, options=["GEO_ENCODING=OGC_BP", "DPI=300"] - ) - tst.testCreateCopy( - check_minmax=0, - check_gt=1, - check_srs=True, - check_checksum_not_null=pdf_checksum_available(), - ) - - -def test_pdf_ogcbp_lcc(poppler_or_pdfium): - wkt = """PROJCS["NAD83 / Utah North", - GEOGCS["NAD83", - DATUM["North_American_Datum_1983", - SPHEROID["GRS 1980",6378137,298.257222101, - AUTHORITY["EPSG","7019"]], - TOWGS84[0,0,0,0,0,0,0]], - PRIMEM["Greenwich",0], - UNIT["degree",0.0174532925199433]], - PROJECTION["Lambert_Conformal_Conic_2SP"], - PARAMETER["standard_parallel_1",41.78333333333333], - PARAMETER["standard_parallel_2",40.71666666666667], - PARAMETER["latitude_of_origin",40.33333333333334], - PARAMETER["central_meridian",-111.5], - PARAMETER["false_easting",500000], - PARAMETER["false_northing",1000000], - UNIT["metre",1]]]""" - - src_ds = gdal.GetDriverByName("GTiff").Create("tmp/temp.tif", 1, 1) - src_ds.SetProjection(wkt) - src_ds.SetGeoTransform([500000, 1, 0, 1000000, 0, -1]) - - with gdal.config_option("GDAL_PDF_OGC_BP_WRITE_WKT", "FALSE"): - out_ds = gdaltest.pdf_drv.CreateCopy("tmp/pdf_ogcbp_lcc.pdf", src_ds) - out_wkt = out_ds.GetProjectionRef() - out_ds = None - - src_ds = None - - gdal.Unlink("tmp/temp.tif") - gdal.GetDriverByName("PDF").Delete("tmp/pdf_ogcbp_lcc.pdf") - - sr1 = osr.SpatialReference(wkt) - sr2 = osr.SpatialReference(out_wkt) - if sr1.IsSame(sr2) == 0: - print(sr2.ExportToPrettyWkt()) - pytest.fail("wrong wkt") - - ############################################################################### # Test no compression @@ -1001,138 +926,6 @@ def test_pdf_update_gcps_iso32000(poppler_or_pdfium): _pdf_update_gcps(poppler_or_pdfium) -@pytest.mark.skipif( - not gdaltest.vrt_has_open_support(), - reason="VRT driver open missing", -) -def test_pdf_update_gcps_ogc_bp(poppler_or_pdfium): - with gdal.config_option("GDAL_PDF_GEO_ENCODING", "OGC_BP"): - _pdf_update_gcps(poppler_or_pdfium) - - -############################################################################### -# Check SetGCPs() but with GCPs that do *not* resolve to a geotransform - - -@pytest.mark.skipif( - not gdaltest.vrt_has_open_support(), - reason="VRT driver open missing", -) -def test_pdf_set_5_gcps_ogc_bp(poppler_or_pdfium): - dpi = 300 - out_filename = "tmp/pdf_set_5_gcps_ogc_bp.pdf" - - src_ds = gdal.Open("data/byte.tif") - src_wkt = src_ds.GetProjectionRef() - src_gt = src_ds.GetGeoTransform() - src_ds = None - - gcp = [ - [2.0, 8.0, 0, 0], - [2.0, 10.0, 0, 0], - [2.0, 18.0, 0, 0], - [16.0, 18.0, 0, 0], - [16.0, 8.0, 0, 0], - ] - - for i, _ in enumerate(gcp): - gcp[i][2] = src_gt[0] + gcp[i][0] * src_gt[1] + gcp[i][1] * src_gt[2] - gcp[i][3] = src_gt[3] + gcp[i][0] * src_gt[4] + gcp[i][1] * src_gt[5] - - # That way, GCPs will not resolve to a geotransform - gcp[1][2] -= 100 - - vrt_txt = """ - - - - - - - - - - data/byte.tif - - 1 - - -""" % ( - src_wkt, - gcp[0][0], - gcp[0][1], - gcp[0][2], - gcp[0][3], - gcp[1][0], - gcp[1][1], - gcp[1][2], - gcp[1][3], - gcp[2][0], - gcp[2][1], - gcp[2][2], - gcp[2][3], - gcp[3][0], - gcp[3][1], - gcp[3][2], - gcp[3][3], - gcp[4][0], - gcp[4][1], - gcp[4][2], - gcp[4][3], - ) - vrt_ds = gdal.Open(vrt_txt) - vrt_gcps = vrt_ds.GetGCPs() - - # Create PDF - ds = gdaltest.pdf_drv.CreateCopy( - out_filename, vrt_ds, options=["GEO_ENCODING=OGC_BP", "DPI=%d" % dpi] - ) - ds = None - - vrt_ds = None - - # Check - ds = gdal.Open(out_filename) - got_gt = ds.GetGeoTransform() - got_wkt = ds.GetProjectionRef() - got_gcp_count = ds.GetGCPCount() - got_gcps = ds.GetGCPs() - got_gcp_wkt = ds.GetGCPProjection() - got_neatline = ds.GetMetadataItem("NEATLINE") - ds = None - - assert got_wkt == "", "did not expect non null GetProjectionRef" - - assert got_gcp_wkt != "", "did not expect null GetGCPProjection" - - expected_gt = [0, 1, 0, 0, 0, 1] - for i in range(6): - assert got_gt[i] == pytest.approx( - expected_gt[i], abs=1e-8 - ), "did not get expected gt" - - assert got_gcp_count == len(gcp), "did not get expected GCP count" - - for i in range(got_gcp_count): - assert ( - got_gcps[i].GCPX == pytest.approx(vrt_gcps[i].GCPX, abs=1e-5) - and got_gcps[i].GCPY == pytest.approx(vrt_gcps[i].GCPY, abs=1e-5) - and got_gcps[i].GCPPixel == pytest.approx(vrt_gcps[i].GCPPixel, abs=1e-5) - and got_gcps[i].GCPLine == pytest.approx(vrt_gcps[i].GCPLine, abs=1e-5) - ), ("did not get expected GCP (%d)" % i) - - got_geom = ogr.CreateGeometryFromWkt(got_neatline) - # Not sure this is really what we want, but without any geotransform, we cannot - # find projected coordinates - expected_geom = ogr.CreateGeometryFromWkt( - "POLYGON ((2 8,2 10,2 18,16 18,16 8,2 8))" - ) - - ogrtest.check_feature_geometry(got_geom, expected_geom) - - gdaltest.pdf_drv.Delete(out_filename) - - ############################################################################### # Check NEATLINE support @@ -1235,10 +1028,6 @@ def test_pdf_set_neatline_iso32000(poppler_or_pdfium): return _pdf_set_neatline(poppler_or_pdfium, "ISO32000") -def test_pdf_set_neatline_ogc_bp(poppler_or_pdfium): - return _pdf_set_neatline(poppler_or_pdfium, "OGC_BP") - - ############################################################################### # Check that we can generate identical file @@ -1276,46 +1065,6 @@ def test_pdf_check_identity_iso32000(poppler_or_pdfium): gdaltest.pdf_drv.Delete(out_filename) -############################################################################### -# Check that we can generate identical file - - -@pytest.mark.skipif( - not gdaltest.vrt_has_open_support(), - reason="VRT driver open missing", -) -def test_pdf_check_identity_ogc_bp(poppler_or_pdfium): - out_filename = "tmp/pdf_check_identity_ogc_bp.pdf" - - src_ds = gdal.Open("data/pdf/test_pdf.vrt") - with gdal.config_option("GDAL_PDF_OGC_BP_WRITE_WKT", "NO"): - out_ds = gdaltest.pdf_drv.CreateCopy( - out_filename, - src_ds, - options=["GEO_ENCODING=OGC_BP", "STREAM_COMPRESS=NONE"], - ) - del out_ds - src_ds = None - - f = open("data/pdf/test_ogc_bp.pdf", "rb") - data_ref = f.read() - f.close() - - f = open("data/pdf/test_ogc_bp_libpng_1_6_40.pdf", "rb") - data_ref2 = f.read() - f.close() - - f = open(out_filename, "rb") - data_got = f.read() - f.close() - - assert ( - data_got == data_ref or data_got == data_ref2 - ), "content does not match reference content" - - gdaltest.pdf_drv.Delete(out_filename) - - ############################################################################### # Check layers support @@ -2093,7 +1842,7 @@ def test_pdf_composition(): 20 10 - + EPSG:4326 POLYGON((1 1,19 1,19 9,1 9,1 1)) diff --git a/doc/source/drivers/raster/pdf.rst b/doc/source/drivers/raster/pdf.rst index 92a53c35bfd3..eb3acf806be2 100644 --- a/doc/source/drivers/raster/pdf.rst +++ b/doc/source/drivers/raster/pdf.rst @@ -18,8 +18,8 @@ on top of the raster layer (see OGR\_\* creation options in the below section). The driver supports reading georeferencing encoded in either of the 2 -current existing ways : according to the OGC encoding best practice, or -according to the Adobe Supplement to ISO 32000. +current existing ways : according to the Adobe Supplement to ISO 32000, +or to a so-called OGC encoding "best practice" (which is not recommended at all). Multipage documents are exposed as subdatasets, one subdataset par page of the document. @@ -211,9 +211,7 @@ PDF documents can be created from other GDAL raster datasets, that have 1 band (graylevel or with color table), 3 bands (RGB) or 4 bands (RGBA). Georeferencing information will be written by default according to the -ISO32000 specification. It is also possible to write it according to the -OGC Best Practice conventions (but limited to a few datum and projection -types). +ISO32000 specification. Note: PDF write support does not require linking to any backend. @@ -381,11 +379,13 @@ The following creation options are available: Margin below image in user units. - .. co:: GEO_ENCODING - :choices: NONE, ISO32000, OGC_BP, BOTH + :choices: NONE, ISO32000 :default: ISO32000 Set the Geo encoding method to use. + Support for deprecated OGC_BP and BOTH options has been removed in GDAL 3.11 + - .. co:: NEATLINE :choices: @@ -536,8 +536,7 @@ mode in order to set or update the following elements : - xml:XMP metadata (with SetMetadata(md, "xml:XMP")) For geotransform or GCPs, the Geo encoding method used by default is -ISO32000. OGC_BP can be selected by setting the GDAL_PDF_GEO_ENCODING -configuration option to OGC_BP. +ISO32000. Updated elements are written at the end of the file, following the incremental update method described in the PDF specification. @@ -809,14 +808,14 @@ See also Specifications : -- `OGC PDF Georegistration Encoding Best Practice Version 2.2 - (08-139r3) `__ - `Adobe Supplement to ISO 32000 `__ - `PDF Reference, version 1.7 `__ - `Acrobat(R) JavaScript Scripting Reference `__ +- `OGC PDF Georegistration Encoding Best Practice Version 2.2 + (08-139r3) `__ (not recommended) Libraries : diff --git a/frmts/pdf/pdfcreatecopy.cpp b/frmts/pdf/pdfcreatecopy.cpp index 3dd8c18627aa..ec0d25c0c622 100644 --- a/frmts/pdf/pdfcreatecopy.cpp +++ b/frmts/pdf/pdfcreatecopy.cpp @@ -325,9 +325,6 @@ void GDALPDFUpdateWriter::UpdateProj(GDALDataset *poSrcDS, double dfDPI, if (EQUAL(pszGEO_ENCODING, "ISO32000") || EQUAL(pszGEO_ENCODING, "BOTH")) nViewportId = WriteSRS_ISO32000(poSrcDS, dfDPI * USER_UNIT_IN_INCH, nullptr, &sMargins, TRUE); - if (EQUAL(pszGEO_ENCODING, "OGC_BP") || EQUAL(pszGEO_ENCODING, "BOTH")) - nLGIDictId = WriteSRS_OGC_BP(poSrcDS, dfDPI * USER_UNIT_IN_INCH, - nullptr, &sMargins); #ifdef invalidate_xref_entry GDALPDFObject *poVP = poPageDict->Get("VP"); @@ -950,630 +947,6 @@ GDALPDFObjectNum GDALPDFBaseWriter::WriteSRS_ISO32000(GDALDataset *poSrcDS, return nViewportId.toBool() ? nViewportId : nMeasureId; } -/************************************************************************/ -/* GDALPDFBuildOGC_BP_Datum() */ -/************************************************************************/ - -static GDALPDFObject *GDALPDFBuildOGC_BP_Datum(const OGRSpatialReference *poSRS) -{ - const OGR_SRSNode *poDatumNode = poSRS->GetAttrNode("DATUM"); - const char *pszDatumDescription = nullptr; - if (poDatumNode && poDatumNode->GetChildCount() > 0) - pszDatumDescription = poDatumNode->GetChild(0)->GetValue(); - - GDALPDFObjectRW *poPDFDatum = nullptr; - - if (pszDatumDescription) - { - double dfSemiMajor = poSRS->GetSemiMajor(); - double dfInvFlattening = poSRS->GetInvFlattening(); - int nEPSGDatum = -1; - const char *pszAuthority = poSRS->GetAuthorityName("DATUM"); - if (pszAuthority != nullptr && EQUAL(pszAuthority, "EPSG")) - nEPSGDatum = atoi(poSRS->GetAuthorityCode("DATUM")); - - if (EQUAL(pszDatumDescription, SRS_DN_WGS84) || nEPSGDatum == 6326) - poPDFDatum = GDALPDFObjectRW::CreateString("WGE"); - else if (EQUAL(pszDatumDescription, SRS_DN_NAD27) || nEPSGDatum == 6267) - poPDFDatum = GDALPDFObjectRW::CreateString("NAS"); - else if (EQUAL(pszDatumDescription, SRS_DN_NAD83) || nEPSGDatum == 6269) - poPDFDatum = GDALPDFObjectRW::CreateString("NAR"); - else if (nEPSGDatum == 6135) - poPDFDatum = GDALPDFObjectRW::CreateString("OHA-M"); - else - { - CPLDebug("PDF", - "Unhandled datum name (%s). Write datum parameters then.", - pszDatumDescription); - - GDALPDFDictionaryRW *poPDFDatumDict = new GDALPDFDictionaryRW(); - poPDFDatum = GDALPDFObjectRW::CreateDictionary(poPDFDatumDict); - - const OGR_SRSNode *poSpheroidNode = poSRS->GetAttrNode("SPHEROID"); - if (poSpheroidNode && poSpheroidNode->GetChildCount() >= 3) - { - poPDFDatumDict->Add("Description", pszDatumDescription); - -#ifdef disabled_because_terrago_toolbar_does_not_like_it - const char *pszEllipsoidCode = NULL; - if (std::abs(dfSemiMajor - 6378249.145) < 0.01 && - std::abs(dfInvFlattening - 293.465) < 0.0001) - { - pszEllipsoidCode = "CD"; /* Clark 1880 */ - } - else if (std::abs(dfSemiMajor - 6378245.0) < 0.01 x && - std::abs(dfInvFlattening - 298.3) < 0.0001) - { - pszEllipsoidCode = "KA"; /* Krassovsky */ - } - else if (std::abs(dfSemiMajor - 6378388.0) < 0.01 && - std::abs(dfInvFlattening - 297.0) < 0.0001) - { - pszEllipsoidCode = "IN"; /* International 1924 */ - } - else if (std::abs(dfSemiMajor - 6378160.0) < 0.01 && - std::abs(dfInvFlattening - 298.25) < 0.0001) - { - pszEllipsoidCode = "AN"; /* Australian */ - } - else if (std::abs(dfSemiMajor - 6377397.155) < 0.01 && - std::abs(dfInvFlattening - 299.1528128) < 0.0001) - { - pszEllipsoidCode = "BR"; /* Bessel 1841 */ - } - else if (std::abs(dfSemiMajor - 6377483.865) < 0.01 && - std::abs(dfInvFlattening - 299.1528128) < 0.0001) - { - pszEllipsoidCode = - "BN"; /* Bessel 1841 (Namibia / Schwarzeck)*/ - } -#if 0 - else if( std::abs(dfSemiMajor-6378160.0) < 0.01 - && std::abs(dfInvFlattening-298.247167427) < 0.0001 ) - { - pszEllipsoidCode = "GRS67"; /* GRS 1967 */ - } -#endif - else if (std::abs(dfSemiMajor - 6378137) < 0.01 && - std::abs(dfInvFlattening - 298.257222101) < 0.000001) - { - pszEllipsoidCode = "RF"; /* GRS 1980 */ - } - else if (std::abs(dfSemiMajor - 6378206.4) < 0.01 && - std::abs(dfInvFlattening - 294.9786982) < 0.0001) - { - pszEllipsoidCode = "CC"; /* Clarke 1866 */ - } - else if (std::abs(dfSemiMajor - 6377340.189) < 0.01 && - std::abs(dfInvFlattening - 299.3249646) < 0.0001) - { - pszEllipsoidCode = "AM"; /* Modified Airy */ - } - else if (std::abs(dfSemiMajor - 6377563.396) < 0.01 && - std::abs(dfInvFlattening - 299.3249646) < 0.0001) - { - pszEllipsoidCode = "AA"; /* Airy */ - } - else if (std::abs(dfSemiMajor - 6378200) < 0.01 && - std::abs(dfInvFlattening - 298.3) < 0.0001) - { - pszEllipsoidCode = "HE"; /* Helmert 1906 */ - } - else if (std::abs(dfSemiMajor - 6378155) < 0.01 && - std::abs(dfInvFlattening - 298.3) < 0.0001) - { - pszEllipsoidCode = "FA"; /* Modified Fischer 1960 */ - } -#if 0 - else if( std::abs(dfSemiMajor-6377298.556) < 0.01 - && std::abs(dfInvFlattening-300.8017) < 0.0001 ) - { - pszEllipsoidCode = "evrstSS"; /* Everest (Sabah & Sarawak) */ - } - else if( std::abs(dfSemiMajor-6378165.0) < 0.01 - && std::abs(dfInvFlattening-298.3) < 0.0001 ) - { - pszEllipsoidCode = "WGS60"; - } - else if( std::abs(dfSemiMajor-6378145.0) < 0.01 - && std::abs(dfInvFlattening-298.25) < 0.0001 ) - { - pszEllipsoidCode = "WGS66"; - } -#endif - else if (std::abs(dfSemiMajor - 6378135.0) < 0.01 && - std::abs(dfInvFlattening - 298.26) < 0.0001) - { - pszEllipsoidCode = "WD"; - } - else if (std::abs(dfSemiMajor - 6378137.0) < 0.01 && - std::abs(dfInvFlattening - 298.257223563) < 0.000001) - { - pszEllipsoidCode = "WE"; - } - - if (pszEllipsoidCode != NULL) - { - poPDFDatumDict->Add("Ellipsoid", pszEllipsoidCode); - } - else -#endif /* disabled_because_terrago_toolbar_does_not_like_it */ - { - const char *pszEllipsoidDescription = - poSpheroidNode->GetChild(0)->GetValue(); - - CPLDebug("PDF", - "Unhandled ellipsoid name (%s). Write ellipsoid " - "parameters then.", - pszEllipsoidDescription); - - poPDFDatumDict->Add( - "Ellipsoid", - &((new GDALPDFDictionaryRW()) - ->Add("Description", pszEllipsoidDescription) - .Add("SemiMajorAxis", dfSemiMajor, TRUE) - .Add("InvFlattening", dfInvFlattening, TRUE))); - } - - const OGR_SRSNode *poTOWGS84 = poSRS->GetAttrNode("TOWGS84"); - if (poTOWGS84 != nullptr && poTOWGS84->GetChildCount() >= 3 && - (poTOWGS84->GetChildCount() < 7 || - (EQUAL(poTOWGS84->GetChild(3)->GetValue(), "") && - EQUAL(poTOWGS84->GetChild(4)->GetValue(), "") && - EQUAL(poTOWGS84->GetChild(5)->GetValue(), "") && - EQUAL(poTOWGS84->GetChild(6)->GetValue(), "")))) - { - poPDFDatumDict->Add( - "ToWGS84", - &((new GDALPDFDictionaryRW()) - ->Add("dx", poTOWGS84->GetChild(0)->GetValue()) - .Add("dy", poTOWGS84->GetChild(1)->GetValue()) - .Add("dz", poTOWGS84->GetChild(2)->GetValue()))); - } - else if (poTOWGS84 != nullptr && - poTOWGS84->GetChildCount() >= 7) - { - poPDFDatumDict->Add( - "ToWGS84", - &((new GDALPDFDictionaryRW()) - ->Add("dx", poTOWGS84->GetChild(0)->GetValue()) - .Add("dy", poTOWGS84->GetChild(1)->GetValue()) - .Add("dz", poTOWGS84->GetChild(2)->GetValue()) - .Add("rx", poTOWGS84->GetChild(3)->GetValue()) - .Add("ry", poTOWGS84->GetChild(4)->GetValue()) - .Add("rz", poTOWGS84->GetChild(5)->GetValue()) - .Add("sf", poTOWGS84->GetChild(6)->GetValue()))); - } - } - } - } - else - { - CPLError(CE_Warning, CPLE_NotSupported, - "No datum name. Defaulting to WGS84."); - } - - if (poPDFDatum == nullptr) - poPDFDatum = GDALPDFObjectRW::CreateString("WGE"); - - return poPDFDatum; -} - -/************************************************************************/ -/* GDALPDFBuildOGC_BP_Projection() */ -/************************************************************************/ - -GDALPDFDictionaryRW *GDALPDFBaseWriter::GDALPDFBuildOGC_BP_Projection( - const OGRSpatialReference *poSRS) -{ - - const char *pszProjectionOGCBP = "GEOGRAPHIC"; - const char *pszProjection = poSRS->GetAttrValue("PROJECTION"); - - GDALPDFDictionaryRW *poProjectionDict = new GDALPDFDictionaryRW(); - poProjectionDict->Add("Type", GDALPDFObjectRW::CreateName("Projection")); - poProjectionDict->Add("Datum", GDALPDFBuildOGC_BP_Datum(poSRS)); - - if (pszProjection == nullptr) - { - if (poSRS->IsGeographic()) - pszProjectionOGCBP = "GEOGRAPHIC"; - else if (poSRS->IsLocal()) - pszProjectionOGCBP = "LOCAL CARTESIAN"; - else - { - CPLError(CE_Warning, CPLE_NotSupported, "Unsupported SRS type"); - delete poProjectionDict; - return nullptr; - } - } - else if (EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR)) - { - int bNorth; - int nZone = poSRS->GetUTMZone(&bNorth); - - if (nZone != 0) - { - pszProjectionOGCBP = "UT"; - poProjectionDict->Add("Hemisphere", (bNorth) ? "N" : "S"); - poProjectionDict->Add("Zone", nZone); - } - else - { - double dfCenterLat = - poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 90.L); - double dfCenterLong = - poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0); - double dfScale = poSRS->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0); - double dfFalseEasting = - poSRS->GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0); - double dfFalseNorthing = - poSRS->GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0); - - /* OGC_BP supports representing numbers as strings for better - * precision */ - /* so use it */ - - pszProjectionOGCBP = "TC"; - poProjectionDict->Add("OriginLatitude", dfCenterLat, TRUE); - poProjectionDict->Add("CentralMeridian", dfCenterLong, TRUE); - poProjectionDict->Add("ScaleFactor", dfScale, TRUE); - poProjectionDict->Add("FalseEasting", dfFalseEasting, TRUE); - poProjectionDict->Add("FalseNorthing", dfFalseNorthing, TRUE); - } - } - else if (EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC)) - { - double dfCenterLat = - poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0); - double dfCenterLong = - poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0); - double dfScale = poSRS->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0); - double dfFalseEasting = - poSRS->GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0); - double dfFalseNorthing = - poSRS->GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0); - - if (fabs(dfCenterLat) == 90.0 && dfCenterLong == 0.0 && - dfScale == 0.994 && dfFalseEasting == 200000.0 && - dfFalseNorthing == 200000.0) - { - pszProjectionOGCBP = "UP"; - poProjectionDict->Add("Hemisphere", (dfCenterLat > 0) ? "N" : "S"); - } - else - { - pszProjectionOGCBP = "PG"; - poProjectionDict->Add("LatitudeTrueScale", dfCenterLat, TRUE); - poProjectionDict->Add("LongitudeDownFromPole", dfCenterLong, TRUE); - poProjectionDict->Add("ScaleFactor", dfScale, TRUE); - poProjectionDict->Add("FalseEasting", dfFalseEasting, TRUE); - poProjectionDict->Add("FalseNorthing", dfFalseNorthing, TRUE); - } - } - - else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP)) - { - double dfStdP1 = - poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0); - double dfStdP2 = - poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0); - double dfCenterLat = - poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0); - double dfCenterLong = - poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0); - double dfFalseEasting = - poSRS->GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0); - double dfFalseNorthing = - poSRS->GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0); - - pszProjectionOGCBP = "LE"; - poProjectionDict->Add("StandardParallelOne", dfStdP1, TRUE); - poProjectionDict->Add("StandardParallelTwo", dfStdP2, TRUE); - poProjectionDict->Add("OriginLatitude", dfCenterLat, TRUE); - poProjectionDict->Add("CentralMeridian", dfCenterLong, TRUE); - poProjectionDict->Add("FalseEasting", dfFalseEasting, TRUE); - poProjectionDict->Add("FalseNorthing", dfFalseNorthing, TRUE); - } - - else if (EQUAL(pszProjection, SRS_PT_MERCATOR_1SP)) - { - double dfCenterLong = - poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0); - double dfCenterLat = - poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0); - double dfScale = poSRS->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0); - double dfFalseEasting = - poSRS->GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0); - double dfFalseNorthing = - poSRS->GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0); - - pszProjectionOGCBP = "MC"; - poProjectionDict->Add("CentralMeridian", dfCenterLong, TRUE); - poProjectionDict->Add("OriginLatitude", dfCenterLat, TRUE); - poProjectionDict->Add("ScaleFactor", dfScale, TRUE); - poProjectionDict->Add("FalseEasting", dfFalseEasting, TRUE); - poProjectionDict->Add("FalseNorthing", dfFalseNorthing, TRUE); - } - -#ifdef not_supported - else if (EQUAL(pszProjection, SRS_PT_MERCATOR_2SP)) - { - double dfStdP1 = - poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0); - double dfCenterLong = - poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0); - double dfFalseEasting = - poSRS->GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0); - double dfFalseNorthing = - poSRS->GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0); - - pszProjectionOGCBP = "MC"; - poProjectionDict->Add("StandardParallelOne", dfStdP1, TRUE); - poProjectionDict->Add("CentralMeridian", dfCenterLong, TRUE); - poProjectionDict->Add("FalseEasting", dfFalseEasting, TRUE); - poProjectionDict->Add("FalseNorthing", dfFalseNorthing, TRUE); - } -#endif - - else - { - CPLError(CE_Warning, CPLE_NotSupported, - "Unhandled projection type (%s) for now", pszProjection); - } - - poProjectionDict->Add("ProjectionType", pszProjectionOGCBP); - - if (poSRS->IsProjected()) - { - const char *pszUnitName = nullptr; - double dfLinearUnits = poSRS->GetLinearUnits(&pszUnitName); - if (dfLinearUnits == 1.0) - poProjectionDict->Add("Units", "M"); - else if (dfLinearUnits == 0.3048) - poProjectionDict->Add("Units", "FT"); - } - - return poProjectionDict; -} - -/************************************************************************/ -/* WriteSRS_OGC_BP() */ -/************************************************************************/ - -GDALPDFObjectNum GDALPDFBaseWriter::WriteSRS_OGC_BP(GDALDataset *poSrcDS, - double dfUserUnit, - const char *pszNEATLINE, - PDFMargins *psMargins) -{ - int nWidth = poSrcDS->GetRasterXSize(); - int nHeight = poSrcDS->GetRasterYSize(); - const char *pszWKT = poSrcDS->GetProjectionRef(); - double adfGeoTransform[6]; - - int bHasGT = (poSrcDS->GetGeoTransform(adfGeoTransform) == CE_None); - int nGCPCount = poSrcDS->GetGCPCount(); - const GDAL_GCP *pasGCPList = - (nGCPCount >= 4) ? poSrcDS->GetGCPs() : nullptr; - if (pasGCPList != nullptr) - pszWKT = poSrcDS->GetGCPProjection(); - - if (!bHasGT && pasGCPList == nullptr) - return GDALPDFObjectNum(); - - if (pszWKT == nullptr || EQUAL(pszWKT, "")) - return GDALPDFObjectNum(); - - if (!bHasGT) - { - if (!GDALGCPsToGeoTransform(nGCPCount, pasGCPList, adfGeoTransform, - FALSE)) - { - CPLDebug("PDF", "Could not compute GT with exact match. Writing " - "Registration then"); - } - else - { - bHasGT = TRUE; - } - } - - OGRSpatialReferenceH hSRS = OSRNewSpatialReference(pszWKT); - if (hSRS == nullptr) - return GDALPDFObjectNum(); - OSRSetAxisMappingStrategy(hSRS, OAMS_TRADITIONAL_GIS_ORDER); - - const OGRSpatialReference *poSRS = OGRSpatialReference::FromHandle(hSRS); - GDALPDFDictionaryRW *poProjectionDict = - GDALPDFBuildOGC_BP_Projection(poSRS); - if (poProjectionDict == nullptr) - { - OSRDestroySpatialReference(hSRS); - return GDALPDFObjectNum(); - } - - GDALPDFArrayRW *poNeatLineArray = nullptr; - - if (pszNEATLINE == nullptr) - pszNEATLINE = poSrcDS->GetMetadataItem("NEATLINE"); - if (bHasGT && pszNEATLINE != nullptr && !EQUAL(pszNEATLINE, "NO") && - pszNEATLINE[0] != '\0') - { - OGRGeometry *poGeom = nullptr; - OGRGeometryFactory::createFromWkt(pszNEATLINE, nullptr, &poGeom); - if (poGeom != nullptr && - wkbFlatten(poGeom->getGeometryType()) == wkbPolygon) - { - OGRLineString *poLS = poGeom->toPolygon()->getExteriorRing(); - double adfGeoTransformInv[6]; - if (poLS != nullptr && poLS->getNumPoints() >= 5 && - GDALInvGeoTransform(adfGeoTransform, adfGeoTransformInv)) - { - poNeatLineArray = new GDALPDFArrayRW(); - - // FIXME : ensure that they are in clockwise order ? - for (int i = 0; i < poLS->getNumPoints() - 1; i++) - { - double X = poLS->getX(i); - double Y = poLS->getY(i); - double x = adfGeoTransformInv[0] + - X * adfGeoTransformInv[1] + - Y * adfGeoTransformInv[2]; - double y = adfGeoTransformInv[3] + - X * adfGeoTransformInv[4] + - Y * adfGeoTransformInv[5]; - poNeatLineArray->Add(x / dfUserUnit + psMargins->nLeft, - TRUE); - poNeatLineArray->Add( - (nHeight - y) / dfUserUnit + psMargins->nBottom, TRUE); - } - } - } - delete poGeom; - } - - if (pszNEATLINE != nullptr && EQUAL(pszNEATLINE, "NO")) - { - // Do nothing - } - else if (pasGCPList && poNeatLineArray == nullptr) - { - if (nGCPCount == 4) - { - int iUL = 0; - int iUR = 0; - int iLR = 0; - int iLL = 0; - GDALPDFFind4Corners(pasGCPList, iUL, iUR, iLR, iLL); - - double adfNL[8]; - adfNL[0] = - pasGCPList[iUL].dfGCPPixel / dfUserUnit + psMargins->nLeft; - adfNL[1] = (nHeight - pasGCPList[iUL].dfGCPLine) / dfUserUnit + - psMargins->nBottom; - adfNL[2] = - pasGCPList[iLL].dfGCPPixel / dfUserUnit + psMargins->nLeft; - adfNL[3] = (nHeight - pasGCPList[iLL].dfGCPLine) / dfUserUnit + - psMargins->nBottom; - adfNL[4] = - pasGCPList[iLR].dfGCPPixel / dfUserUnit + psMargins->nLeft; - adfNL[5] = (nHeight - pasGCPList[iLR].dfGCPLine) / dfUserUnit + - psMargins->nBottom; - adfNL[6] = - pasGCPList[iUR].dfGCPPixel / dfUserUnit + psMargins->nLeft; - adfNL[7] = (nHeight - pasGCPList[iUR].dfGCPLine) / dfUserUnit + - psMargins->nBottom; - - poNeatLineArray = new GDALPDFArrayRW(); - poNeatLineArray->Add(adfNL, 8, TRUE); - } - else - { - poNeatLineArray = new GDALPDFArrayRW(); - - // FIXME : ensure that they are in clockwise order ? - int i; - for (i = 0; i < nGCPCount; i++) - { - poNeatLineArray->Add(pasGCPList[i].dfGCPPixel / dfUserUnit + - psMargins->nLeft, - TRUE); - poNeatLineArray->Add((nHeight - pasGCPList[i].dfGCPLine) / - dfUserUnit + - psMargins->nBottom, - TRUE); - } - } - } - else if (poNeatLineArray == nullptr) - { - poNeatLineArray = new GDALPDFArrayRW(); - - poNeatLineArray->Add(0 / dfUserUnit + psMargins->nLeft, TRUE); - poNeatLineArray->Add((nHeight - 0) / dfUserUnit + psMargins->nBottom, - TRUE); - - poNeatLineArray->Add(0 / dfUserUnit + psMargins->nLeft, TRUE); - poNeatLineArray->Add( - (/*nHeight -nHeight*/ 0) / dfUserUnit + psMargins->nBottom, TRUE); - - poNeatLineArray->Add(nWidth / dfUserUnit + psMargins->nLeft, TRUE); - poNeatLineArray->Add( - (/*nHeight -nHeight*/ 0) / dfUserUnit + psMargins->nBottom, TRUE); - - poNeatLineArray->Add(nWidth / dfUserUnit + psMargins->nLeft, TRUE); - poNeatLineArray->Add((nHeight - 0) / dfUserUnit + psMargins->nBottom, - TRUE); - } - - auto nLGIDictId = AllocNewObject(); - StartObj(nLGIDictId); - GDALPDFDictionaryRW oLGIDict; - oLGIDict.Add("Type", GDALPDFObjectRW::CreateName("LGIDict")) - .Add("Version", "2.1"); - if (bHasGT) - { - double adfCTM[6]; - double dfX1 = psMargins->nLeft; - double dfY2 = nHeight / dfUserUnit + psMargins->nBottom; - - adfCTM[0] = adfGeoTransform[1] * dfUserUnit; - adfCTM[1] = adfGeoTransform[2] * dfUserUnit; - adfCTM[2] = -adfGeoTransform[4] * dfUserUnit; - adfCTM[3] = -adfGeoTransform[5] * dfUserUnit; - adfCTM[4] = adfGeoTransform[0] - (adfCTM[0] * dfX1 + adfCTM[2] * dfY2); - adfCTM[5] = adfGeoTransform[3] - (adfCTM[1] * dfX1 + adfCTM[3] * dfY2); - - oLGIDict.Add("CTM", &((new GDALPDFArrayRW())->Add(adfCTM, 6, TRUE))); - } - else - { - GDALPDFArrayRW *poRegistrationArray = new GDALPDFArrayRW(); - int i; - for (i = 0; i < nGCPCount; i++) - { - GDALPDFArrayRW *poPTArray = new GDALPDFArrayRW(); - poPTArray->Add( - pasGCPList[i].dfGCPPixel / dfUserUnit + psMargins->nLeft, TRUE); - poPTArray->Add((nHeight - pasGCPList[i].dfGCPLine) / dfUserUnit + - psMargins->nBottom, - TRUE); - poPTArray->Add(pasGCPList[i].dfGCPX, TRUE); - poPTArray->Add(pasGCPList[i].dfGCPY, TRUE); - poRegistrationArray->Add(poPTArray); - } - oLGIDict.Add("Registration", poRegistrationArray); - } - if (poNeatLineArray) - { - oLGIDict.Add("Neatline", poNeatLineArray); - } - - const OGR_SRSNode *poNode = poSRS->GetRoot(); - if (poNode != nullptr) - poNode = poNode->GetChild(0); - const char *pszDescription = nullptr; - if (poNode != nullptr) - pszDescription = poNode->GetValue(); - if (pszDescription) - { - oLGIDict.Add("Description", pszDescription); - } - - oLGIDict.Add("Projection", poProjectionDict); - - /* GDAL extension */ - if (CPLTestBool(CPLGetConfigOption("GDAL_PDF_OGC_BP_WRITE_WKT", "TRUE"))) - poProjectionDict->Add("WKT", pszWKT); - - VSIFPrintfL(m_fp, "%s\n", oLGIDict.Serialize().c_str()); - EndObj(); - - OSRDestroySpatialReference(hSRS); - - return nLGIDictId; -} - /************************************************************************/ /* GDALPDFGetValueFromDSOrOption() */ /************************************************************************/ @@ -1758,19 +1131,12 @@ bool GDALPDFWriter::StartPage(GDALDataset *poClippingDS, double dfDPI, const bool bISO32000 = EQUAL(pszGEO_ENCODING, "ISO32000") || EQUAL(pszGEO_ENCODING, "BOTH"); - const bool bOGC_BP = - EQUAL(pszGEO_ENCODING, "OGC_BP") || EQUAL(pszGEO_ENCODING, "BOTH"); GDALPDFObjectNum nViewportId; if (bISO32000) nViewportId = WriteSRS_ISO32000(poClippingDS, dfUserUnit, pszNEATLINE, psMargins, TRUE); - GDALPDFObjectNum nLGIDictId; - if (bOGC_BP) - nLGIDictId = - WriteSRS_OGC_BP(poClippingDS, dfUserUnit, pszNEATLINE, psMargins); - StartObj(nPageId); GDALPDFDictionaryRW oDictPage; oDictPage.Add("Type", GDALPDFObjectRW::CreateName("Page")) @@ -1799,10 +1165,6 @@ bool GDALPDFWriter::StartPage(GDALDataset *poClippingDS, double dfDPI, { oDictPage.Add("VP", &((new GDALPDFArrayRW())->Add(nViewportId, 0))); } - if (nLGIDictId.toBool()) - { - oDictPage.Add("LGIDict", nLGIDictId, 0); - } #ifndef HACK_TO_GENERATE_OCMD if (bHasOGRData) @@ -4912,6 +4274,20 @@ GDALDataset *GDALPDFCreateCopy(const char *pszFilename, GDALDataset *poSrcDS, const char *pszGEO_ENCODING = CSLFetchNameValueDef(papszOptions, "GEO_ENCODING", "ISO32000"); + if (EQUAL(pszGEO_ENCODING, "OGC_BP")) + { + CPLError(CE_Failure, CPLE_NotSupported, + "GEO_ENCODING=OGC_BP is no longer supported. Switch to using " + "ISO32000"); + return nullptr; + } + else if (EQUAL(pszGEO_ENCODING, "BOTH")) + { + CPLError(CE_Warning, CPLE_NotSupported, + "GEO_ENCODING=BOTH is no longer strictly supported. This now " + "fallbacks to ISO32000"); + pszGEO_ENCODING = "ISO32000"; + } const char *pszXMP = CSLFetchNameValue(papszOptions, "XMP"); diff --git a/frmts/pdf/pdfcreatecopy.h b/frmts/pdf/pdfcreatecopy.h index dd657f88c73d..07f367c5f30b 100644 --- a/frmts/pdf/pdfcreatecopy.h +++ b/frmts/pdf/pdfcreatecopy.h @@ -225,11 +225,6 @@ class GDALPDFBaseWriter const char *pszNEATLINE, PDFMargins *psMargins, int bWriteViewport); - GDALPDFObjectNum WriteSRS_OGC_BP(GDALDataset *poSrcDS, double dfUserUnit, - const char *pszNEATLINE, - PDFMargins *psMargins); - static GDALPDFDictionaryRW * - GDALPDFBuildOGC_BP_Projection(const OGRSpatialReference *poSRS); GDALPDFObjectNum WriteOCG(const char *pszLayerName, diff --git a/frmts/pdf/pdfcreatefromcomposition.cpp b/frmts/pdf/pdfcreatefromcomposition.cpp index f245ac5cb84c..84f25903ceb3 100644 --- a/frmts/pdf/pdfcreatefromcomposition.cpp +++ b/frmts/pdf/pdfcreatefromcomposition.cpp @@ -625,7 +625,7 @@ bool GDALPDFComposerWriter::CreateOutline(const CPLXMLNode *psNode) bool GDALPDFComposerWriter::GenerateGeoreferencing( const CPLXMLNode *psGeoreferencing, double dfWidthInUserUnit, double dfHeightInUserUnit, GDALPDFObjectNum &nViewportId, - GDALPDFObjectNum &nLGIDictId, Georeferencing &georeferencing) + Georeferencing &georeferencing) { double bboxX1 = 0; double bboxY1 = 0; @@ -738,13 +738,10 @@ bool GDALPDFComposerWriter::GenerateGeoreferencing( if (CPLTestBool( CPLGetXMLValue(psGeoreferencing, "OGCBestPracticeFormat", "false"))) { - nLGIDictId = GenerateOGC_BP_Georeferencing( - OGRSpatialReference::ToHandle(poSRS.get()), bboxX1, bboxY1, bboxX2, - bboxY2, aGCPs, aBoundingPolygon); - if (!nLGIDictId.toBool()) - { - return false; - } + CPLError(CE_Failure, CPLE_NotSupported, + "OGCBestPracticeFormat no longer supported. Use " + "ISO32000ExtensionFormat"); + return false; } const char *pszId = CPLGetXMLValue(psGeoreferencing, "id", nullptr); @@ -910,76 +907,6 @@ GDALPDFObjectNum GDALPDFComposerWriter::GenerateISO32000_Georeferencing( return nViewportId; } -/************************************************************************/ -/* GenerateOGC_BP_Georeferencing() */ -/************************************************************************/ - -GDALPDFObjectNum GDALPDFComposerWriter::GenerateOGC_BP_Georeferencing( - OGRSpatialReferenceH hSRS, double bboxX1, double bboxY1, double bboxX2, - double bboxY2, const std::vector &aGCPs, - const std::vector &aBoundingPolygon) -{ - const OGRSpatialReference *poSRS = OGRSpatialReference::FromHandle(hSRS); - GDALPDFDictionaryRW *poProjectionDict = - GDALPDFBuildOGC_BP_Projection(poSRS); - if (poProjectionDict == nullptr) - { - OSRDestroySpatialReference(hSRS); - return GDALPDFObjectNum(); - } - - GDALPDFArrayRW *poNeatLineArray = new GDALPDFArrayRW(); - if (!aBoundingPolygon.empty()) - { - for (const auto &xy : aBoundingPolygon) - { - poNeatLineArray->Add(xy.x).Add(xy.y); - } - } - else - { - poNeatLineArray->Add(bboxX1).Add(bboxY1).Add(bboxX2).Add(bboxY2); - } - - GDALPDFArrayRW *poRegistration = new GDALPDFArrayRW(); - - for (const auto &gcp : aGCPs) - { - GDALPDFArrayRW *poGCP = new GDALPDFArrayRW(); - poGCP->Add(gcp.Pixel(), TRUE) - .Add(gcp.Line(), TRUE) - .Add(gcp.X(), TRUE) - .Add(gcp.Y(), TRUE); - poRegistration->Add(poGCP); - } - - auto nLGIDictId = AllocNewObject(); - StartObj(nLGIDictId); - GDALPDFDictionaryRW oLGIDict; - oLGIDict.Add("Type", GDALPDFObjectRW::CreateName("LGIDict")) - .Add("Version", "2.1") - .Add("Neatline", poNeatLineArray); - - oLGIDict.Add("Registration", poRegistration); - - /* GDAL extension */ - if (CPLTestBool(CPLGetConfigOption("GDAL_PDF_OGC_BP_WRITE_WKT", "TRUE"))) - { - char *pszWKT = nullptr; - OSRExportToWkt(hSRS, &pszWKT); - if (pszWKT) - poProjectionDict->Add("WKT", pszWKT); - CPLFree(pszWKT); - } - - oLGIDict.Add("Projection", poProjectionDict); - - VSIFPrintfL(m_fp, "%s\n", oLGIDict.Serialize().c_str()); - EndObj(); - - return nLGIDictId; -} - /************************************************************************/ /* GeneratePage() */ /************************************************************************/ @@ -1000,7 +927,6 @@ bool GDALPDFComposerWriter::GeneratePage(const CPLXMLNode *psPage) USER_UNIT_IN_INCH; std::vector anViewportIds; - std::vector anLGIDictIds; PageContext oPageContext; for (const auto *psIter = psPage->psChild; psIter; psIter = psIter->psNext) @@ -1009,18 +935,15 @@ bool GDALPDFComposerWriter::GeneratePage(const CPLXMLNode *psPage) strcmp(psIter->pszValue, "Georeferencing") == 0) { GDALPDFObjectNum nViewportId; - GDALPDFObjectNum nLGIDictId; Georeferencing georeferencing; if (!GenerateGeoreferencing(psIter, dfWidthInUserUnit, dfHeightInUserUnit, nViewportId, - nLGIDictId, georeferencing)) + georeferencing)) { return false; } if (nViewportId.toBool()) anViewportIds.emplace_back(nViewportId); - if (nLGIDictId.toBool()) - anLGIDictIds.emplace_back(nLGIDictId); if (!georeferencing.m_osID.empty()) { oPageContext.m_oMapGeoreferencedId[georeferencing.m_osID] = @@ -1129,18 +1052,6 @@ bool GDALPDFComposerWriter::GeneratePage(const CPLXMLNode *psPage) oDictPage.Add("VP", poViewports); } - if (anLGIDictIds.size() == 1) - { - oDictPage.Add("LGIDict", anLGIDictIds[0], 0); - } - else if (!anLGIDictIds.empty()) - { - auto poLGIDict = new GDALPDFArrayRW(); - for (const auto &id : anLGIDictIds) - poLGIDict->Add(id, 0); - oDictPage.Add("LGIDict", poLGIDict); - } - if (nStructParentsIdx >= 0) { oDictPage.Add("StructParents", nStructParentsIdx); diff --git a/frmts/pdf/pdfcreatefromcomposition.h b/frmts/pdf/pdfcreatefromcomposition.h index 4c7f72a45f78..4a1ae244b2fc 100644 --- a/frmts/pdf/pdfcreatefromcomposition.h +++ b/frmts/pdf/pdfcreatefromcomposition.h @@ -143,7 +143,6 @@ class GDALPDFComposerWriter final : public GDALPDFBaseWriter double dfWidthInUserUnit, double dfHeightInUserUnit, GDALPDFObjectNum &nViewportId, - GDALPDFObjectNum &nLGIDictId, Georeferencing &georeferencing); GDALPDFObjectNum GenerateISO32000_Georeferencing( @@ -151,12 +150,6 @@ class GDALPDFComposerWriter final : public GDALPDFBaseWriter double bboxY2, const std::vector &aGCPs, const std::vector &aBoundingPolygon); - GDALPDFObjectNum - GenerateOGC_BP_Georeferencing(OGRSpatialReferenceH hSRS, double bboxX1, - double bboxY1, double bboxX2, double bboxY2, - const std::vector &aGCPs, - const std::vector &aBoundingPolygon); - bool ExploreContent(const CPLXMLNode *psNode, PageContext &oPageContext); bool WriteRaster(const CPLXMLNode *psNode, PageContext &oPageContext); bool WriteVector(const CPLXMLNode *psNode, PageContext &oPageContext); diff --git a/frmts/pdf/pdfdrivercore.cpp b/frmts/pdf/pdfdrivercore.cpp index 1ad744359c62..2cf400c959e3 100644 --- a/frmts/pdf/pdfdrivercore.cpp +++ b/frmts/pdf/pdfdrivercore.cpp @@ -150,8 +150,6 @@ void PDFDriverSetCommonMetadata(GDALDriver *poDriver) "description='Format of geo-encoding' default='ISO32000'>\n" " NONE\n" " ISO32000\n" - " OGC_BP\n" - " BOTH\n" " \n" "