|
| 1 | +/****************************************************************************** |
| 2 | + * |
| 3 | + * Project: GDAL |
| 4 | + * Purpose: Common code for gdalalg_raster_clip and gdalalg_vector_clip |
| 5 | + * Author: Even Rouault <even dot rouault at spatialys.com> |
| 6 | + * |
| 7 | + ****************************************************************************** |
| 8 | + * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com> |
| 9 | + * |
| 10 | + * SPDX-License-Identifier: MIT |
| 11 | + ****************************************************************************/ |
| 12 | + |
| 13 | +#include "gdalalg_clip_common.h" |
| 14 | + |
| 15 | +#include "ogrsf_frmts.h" |
| 16 | + |
| 17 | +//! @cond Doxygen_Suppress |
| 18 | + |
| 19 | +#ifndef _ |
| 20 | +#define _(x) (x) |
| 21 | +#endif |
| 22 | + |
| 23 | +/************************************************************************/ |
| 24 | +/* ~GDALClipCommon() */ |
| 25 | +/************************************************************************/ |
| 26 | + |
| 27 | +GDALClipCommon::~GDALClipCommon() = default; |
| 28 | + |
| 29 | +/************************************************************************/ |
| 30 | +/* LoadGeometry() */ |
| 31 | +/************************************************************************/ |
| 32 | + |
| 33 | +std::pair<std::unique_ptr<OGRGeometry>, std::string> |
| 34 | +GDALClipCommon::LoadGeometry() |
| 35 | +{ |
| 36 | + auto poDS = m_likeDataset.GetDatasetRef(); |
| 37 | + OGRLayer *poLyr = nullptr; |
| 38 | + if (!m_likeSQL.empty()) |
| 39 | + poLyr = poDS->ExecuteSQL(m_likeSQL.c_str(), nullptr, nullptr); |
| 40 | + else if (!m_likeLayer.empty()) |
| 41 | + poLyr = poDS->GetLayerByName(m_likeLayer.c_str()); |
| 42 | + else |
| 43 | + poLyr = poDS->GetLayer(0); |
| 44 | + |
| 45 | + if (poLyr == nullptr) |
| 46 | + { |
| 47 | + return {nullptr, |
| 48 | + "Failed to identify source layer from clipping dataset."}; |
| 49 | + } |
| 50 | + |
| 51 | + if (!m_likeWhere.empty()) |
| 52 | + poLyr->SetAttributeFilter(m_likeWhere.c_str()); |
| 53 | + |
| 54 | + OGRGeometryCollection oGC; |
| 55 | + oGC.assignSpatialReference(poLyr->GetSpatialRef()); |
| 56 | + |
| 57 | + for (auto &poFeat : poLyr) |
| 58 | + { |
| 59 | + auto poSrcGeom = std::unique_ptr<OGRGeometry>(poFeat->StealGeometry()); |
| 60 | + if (poSrcGeom) |
| 61 | + { |
| 62 | + // Only take into account areal geometries. |
| 63 | + if (poSrcGeom->getDimension() == 2) |
| 64 | + { |
| 65 | + if (!poSrcGeom->IsValid()) |
| 66 | + { |
| 67 | + return { |
| 68 | + nullptr, |
| 69 | + CPLSPrintf("Geometry of feature " CPL_FRMT_GIB " of %s " |
| 70 | + "is invalid. You may be able to correct it " |
| 71 | + "with 'gdal vector geom make-valid'.", |
| 72 | + poFeat->GetFID(), poDS->GetDescription())}; |
| 73 | + } |
| 74 | + else |
| 75 | + { |
| 76 | + oGC.addGeometry(std::move(poSrcGeom)); |
| 77 | + } |
| 78 | + } |
| 79 | + else |
| 80 | + { |
| 81 | + CPLErrorOnce(CE_Warning, CPLE_AppDefined, |
| 82 | + "Non-polygonal geometry encountered in clipping " |
| 83 | + "dataset will be ignored."); |
| 84 | + } |
| 85 | + } |
| 86 | + } |
| 87 | + |
| 88 | + if (!m_likeSQL.empty()) |
| 89 | + poDS->ReleaseResultSet(poLyr); |
| 90 | + |
| 91 | + if (oGC.IsEmpty()) |
| 92 | + { |
| 93 | + return {nullptr, "No clipping geometry found"}; |
| 94 | + } |
| 95 | + |
| 96 | + return {std::unique_ptr<OGRGeometry>(oGC.UnaryUnion()), std::string()}; |
| 97 | +} |
| 98 | + |
| 99 | +/************************************************************************/ |
| 100 | +/* GetClipGeometry() */ |
| 101 | +/************************************************************************/ |
| 102 | + |
| 103 | +std::pair<std::unique_ptr<OGRGeometry>, std::string> |
| 104 | +GDALClipCommon::GetClipGeometry() |
| 105 | +{ |
| 106 | + |
| 107 | + std::unique_ptr<OGRGeometry> poClipGeom; |
| 108 | + |
| 109 | + if (!m_bbox.empty()) |
| 110 | + { |
| 111 | + poClipGeom = std::make_unique<OGRPolygon>(m_bbox[0], m_bbox[1], |
| 112 | + m_bbox[2], m_bbox[3]); |
| 113 | + |
| 114 | + if (!m_bboxCrs.empty()) |
| 115 | + { |
| 116 | + auto poSRS = new OGRSpatialReference(); |
| 117 | + poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); |
| 118 | + CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput(m_bboxCrs.c_str())); |
| 119 | + poClipGeom->assignSpatialReference(poSRS); |
| 120 | + poSRS->Release(); |
| 121 | + } |
| 122 | + } |
| 123 | + else if (!m_geometry.empty()) |
| 124 | + { |
| 125 | + { |
| 126 | + CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler); |
| 127 | + auto [poGeom, eErr] = |
| 128 | + OGRGeometryFactory::createFromWkt(m_geometry.c_str()); |
| 129 | + if (eErr == OGRERR_NONE) |
| 130 | + { |
| 131 | + poClipGeom = std::move(poGeom); |
| 132 | + } |
| 133 | + else |
| 134 | + { |
| 135 | + poClipGeom.reset( |
| 136 | + OGRGeometryFactory::createFromGeoJson(m_geometry.c_str())); |
| 137 | + if (poClipGeom && poClipGeom->getSpatialReference() == nullptr) |
| 138 | + { |
| 139 | + auto poSRS = new OGRSpatialReference(); |
| 140 | + poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); |
| 141 | + CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput("WGS84")); |
| 142 | + poClipGeom->assignSpatialReference(poSRS); |
| 143 | + poSRS->Release(); |
| 144 | + } |
| 145 | + } |
| 146 | + } |
| 147 | + if (!poClipGeom) |
| 148 | + { |
| 149 | + return { |
| 150 | + nullptr, |
| 151 | + "Clipping geometry is neither a valid WKT or GeoJSON geometry"}; |
| 152 | + } |
| 153 | + |
| 154 | + if (!m_geometryCrs.empty()) |
| 155 | + { |
| 156 | + auto poSRS = new OGRSpatialReference(); |
| 157 | + poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); |
| 158 | + // Validity of CRS already checked by GDALAlgorithm |
| 159 | + CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput(m_geometryCrs.c_str())); |
| 160 | + poClipGeom->assignSpatialReference(poSRS); |
| 161 | + poSRS->Release(); |
| 162 | + } |
| 163 | + } |
| 164 | + else if (auto poLikeDS = m_likeDataset.GetDatasetRef()) |
| 165 | + { |
| 166 | + if (poLikeDS->GetLayerCount() > 1 && m_likeLayer.empty() && |
| 167 | + m_likeSQL.empty()) |
| 168 | + { |
| 169 | + return { |
| 170 | + nullptr, |
| 171 | + "Only single layer dataset can be specified with --like when " |
| 172 | + "neither --like-layer or --like-sql have been specified"}; |
| 173 | + } |
| 174 | + else if (poLikeDS->GetLayerCount() > 0) |
| 175 | + { |
| 176 | + std::string errMsg; |
| 177 | + std::tie(poClipGeom, errMsg) = LoadGeometry(); |
| 178 | + if (!poClipGeom) |
| 179 | + return {nullptr, errMsg}; |
| 180 | + } |
| 181 | + else if (poLikeDS->GetRasterCount() > 0) |
| 182 | + { |
| 183 | + double adfGT[6]; |
| 184 | + if (poLikeDS->GetGeoTransform(adfGT) != CE_None) |
| 185 | + { |
| 186 | + return { |
| 187 | + nullptr, |
| 188 | + CPLSPrintf( |
| 189 | + "Dataset '%s' has no geotransform matrix. Its bounds " |
| 190 | + "cannot be established.", |
| 191 | + poLikeDS->GetDescription())}; |
| 192 | + } |
| 193 | + auto poLikeSRS = poLikeDS->GetSpatialRef(); |
| 194 | + const double dfTLX = adfGT[0]; |
| 195 | + const double dfTLY = adfGT[3]; |
| 196 | + |
| 197 | + double dfTRX = 0; |
| 198 | + double dfTRY = 0; |
| 199 | + GDALApplyGeoTransform(adfGT, poLikeDS->GetRasterXSize(), 0, &dfTRX, |
| 200 | + &dfTRY); |
| 201 | + |
| 202 | + double dfBLX = 0; |
| 203 | + double dfBLY = 0; |
| 204 | + GDALApplyGeoTransform(adfGT, 0, poLikeDS->GetRasterYSize(), &dfBLX, |
| 205 | + &dfBLY); |
| 206 | + |
| 207 | + double dfBRX = 0; |
| 208 | + double dfBRY = 0; |
| 209 | + GDALApplyGeoTransform(adfGT, poLikeDS->GetRasterXSize(), |
| 210 | + poLikeDS->GetRasterYSize(), &dfBRX, &dfBRY); |
| 211 | + |
| 212 | + auto poPoly = std::make_unique<OGRPolygon>(); |
| 213 | + auto poLR = std::make_unique<OGRLinearRing>(); |
| 214 | + poLR->addPoint(dfTLX, dfTLY); |
| 215 | + poLR->addPoint(dfTRX, dfTRY); |
| 216 | + poLR->addPoint(dfBRX, dfBRY); |
| 217 | + poLR->addPoint(dfBLX, dfBLY); |
| 218 | + poLR->addPoint(dfTLX, dfTLY); |
| 219 | + poPoly->addRingDirectly(poLR.release()); |
| 220 | + poPoly->assignSpatialReference(poLikeSRS); |
| 221 | + poClipGeom = std::move(poPoly); |
| 222 | + } |
| 223 | + else |
| 224 | + { |
| 225 | + return {nullptr, "Cannot get extent from clip dataset"}; |
| 226 | + } |
| 227 | + } |
| 228 | + else |
| 229 | + { |
| 230 | + return {nullptr, "--bbox, --geometry or --like must be specified"}; |
| 231 | + } |
| 232 | + |
| 233 | + return {std::move(poClipGeom), std::string()}; |
| 234 | +} |
| 235 | + |
| 236 | +//! @endcond |
0 commit comments