การครอบตัดไฟล์แรสเตอร์ด้วย GDAL warp C ++ API

Aug 20 2020

ฉันต้องการครอบตัดไฟล์แรสเตอร์ DEM ด้วย shapefile Shapefile นี้มีหลายรูปหลายเหลี่ยม ฉันต้องการคำนวณตัวแปรหลายตัวหลังจากการดำเนินการครอบตัด

อัลกอริทึมเขียนด้วย C ++ ดังนั้นฉันต้องการเรียก GDAL C ++ API ฉันไม่จำเป็นต้องเขียนไฟล์วิปริตลงในดิสก์ของฉันและฉันตั้งใจจะใช้ VRT หรือ MEM เมื่อผ่านการทดสอบแล้ว

ฉันได้รับพิกัดของแต่ละรูปหลายเหลี่ยมแล้วและฉันส่งผ่านเป็นรูปหลายเหลี่ยมโดยใช้อาร์กิวเมนต์คัตไลน์ แต่การดำเนินการยังไม่มีผล

            OGRLinearRing poExteriorRing;
            for (i = 0; i < nPt; i++) 
              {
                dX = (*iIterator).vVertex.at(i).dX;
                dY = (*iIterator).vVertex.at(i).dY;
                double dDummy1 = (dX - dX_origin) / dResolution_elevation;
                long lColumn_index = long(round(dDummy1));
                double dDummy2 = (dY_origin - dY) / dResolution_elevation;
                long lRow_index = long(round(dDummy2));     
                poExteriorRing.addPoint(lColumn_index, lRow_index);
              }
            poExteriorRing.closeRings();
            OGRPolygon polygon;
            polygon.addRing(&poExteriorRing);

            GDALDataset  *poVRTDS;
            poDS_elevation->GetGeoTransform(adfGeoTransform);
            poVRTDS = poDriver->CreateCopy( "./test.tiff", poDS_elevation, FALSE, NULL, NULL, NULL );

            GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
            ...
            psWarpOptions->hCutline = &polygon;                
            GDALWarpOperation oOperation;
            oOperation.Initialize( psWarpOptions );
            oOperation.ChunkAndWarpImage( 0, 0,
                            GDALGetRasterXSize( poDS_elevation ),
                            GDALGetRasterYSize( poDS_elevation ) );
            GDALClose( poVRTDS );

ฉันคิดว่าฉันต้องสร้างจุดยอดรูปหลายเหลี่ยมให้กับดัชนีภาพตาม: CUTLINE: สิ่งนี้อาจมีรูปทรงเรขาคณิต WKT สำหรับเส้นตัด จะถูกแปลงเป็นรูปเรขาคณิตโดย GDALWarpOperation :: Initialize () และกำหนดให้กับฟิลด์ GDALWarpOptions hCutline พิกัดต้องแสดงในพิกัดพิกเซล / เส้นต้นทาง หมายเหตุ: สิ่งนี้แตกต่างจากสมมติฐานที่สร้างขึ้นสำหรับตัวเลือก -cutline ของยูทิลิตี้ gdalwarp!https://gdal.org/api/gdalwarp_cpp.html

คำตอบ

1 Chang Aug 24 2020 at 04:20

หลังจากพยายามหลายครั้งฉันสามารถแก้ไขปัญหาหลังจากบล็อกโพสต์บางส่วนได้: http://www.jeepxie.net/article/293437.html

ส่วนหลักคือเราต้องแปลงพิกัด:

eError = TransformCutlineToSource( hWrkSrcDS, hCutline,
                                  &(psWO->papszWarpOptions),
                                  psOptions->papszTO );

และเราไม่ควรใช้ทั้ง hCutline และ cutline ในเวลาเดียวกัน