Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

GDALPolygonize doesn't vectorize all pixels C++ #8001

Closed
sghofrany opened this issue Jun 23, 2023 · 4 comments
Closed

GDALPolygonize doesn't vectorize all pixels C++ #8001

sghofrany opened this issue Jun 23, 2023 · 4 comments

Comments

@sghofrany
Copy link

Expected behavior and actual behavior.

I expect the GDALPolygonize function to take in a raster file and return a vector file that covers all of the pixels in the input.

Steps to reproduce the problem.

Here is the high level of what I am doing:

  1. Read in raster file.
  2. Create a vector layer to pass into GDALPolygonize.
  3. Call GDALPolygonize with raster band from input and the vector layer created.

Here is the code I am running:

    
   string strPath = "chunked_test_chunk_16.tif";

    auto *dataset = (GDALDataset*) GDALOpen(strPath.c_str(), GA_ReadOnly);

    if(dataset == nullptr) {
        return;
    }

    const char *pszDriverName{"GPKG"};
    GDALDriver *poDriver{GetGDALDriverManager()->GetDriverByName(pszDriverName)};

    GDALDataset *polygonized = poDriver->Create((string(filePath.stem()) + ".gpkg").c_str(), 0, 0, 0, GDT_Unknown, nullptr);

    char *cWkt = nullptr;
    dataset->GetSpatialRef()->exportToWkt(&cWkt);

    OGRSpatialReference spatialReference;
    spatialReference.importFromWkt(cWkt);

    OGRwkbGeometryType geomType = wkbPolygon;
    OGRLayer *tmpLayer = polygonized->CreateLayer("polygon", &spatialReference, geomType);

    if(tmpLayer == nullptr) {
        throw std::invalid_argument("Could not create layer");
    }

    OGRFieldDefn field("CellValue", OFTInteger);
    tmpLayer->CreateField(&field);

    GDALRasterBand *rasterBand = dataset->GetRasterBand(1);

    GDALPolygonize(rasterBand, rasterBand, tmpLayer, 0, nullptr, nullptr, nullptr);
    tmpLayer->SyncToDisk();
    
    GDALClose(dataset);
    GDALClose(polygonized);

I have attached the raster file I used and the output file I get from running the above code. GitHub doesn't allow me to upload .tif files, so I have added it in a zip file.

The picture below is what I see in QGIS. The pink color is the polygons and the black is the raster representation. As you can see there the black pixels were not polygonized.

image

Operating system

WSL2
Ubuntu 20.04.2 LTS

GDAL version and provenance

GDAL 3.5.3, released 2022/10/21
chunked_test_chunk_16.zip

@rouault
Copy link
Member

rouault commented Jun 24, 2023

The reason is that you pass the rasterBand itself as the value for the maskBand argument of GDALPolygonize(), which cause pixels at 0 to be considered invalid.
Try rather

GDALPolygonize(rasterBand, rasterBand->GetMaskBand(), tmpLayer, 0, nullptr, nullptr, nullptr);

@rouault rouault closed this as completed Jun 24, 2023
@sghofrany
Copy link
Author

The reason is that you pass the rasterBand itself as the value for the maskBand argument of GDALPolygonize(), which cause pixels at 0 to be considered invalid. Try rather

GDALPolygonize(rasterBand, rasterBand->GetMaskBand(), tmpLayer, 0, nullptr, nullptr, nullptr);

I thought that passing the same raster band again would act as the mask band?

@rouault
Copy link
Member

rouault commented Jun 24, 2023

I thought that passing the same raster band again would act as the mask band?

nothing in the docs suggest so: https://gdal.org/api/gdal_alg.html#_CPPv414GDALPolygonize15GDALRasterBandH15GDALRasterBandH9OGRLayerHiPPc16GDALProgressFuncPv

@sghofrany
Copy link
Author

I thought that passing the same raster band again would act as the mask band?

nothing in the docs suggest so: https://gdal.org/api/gdal_alg.html#_CPPv414GDALPolygonize15GDALRasterBandH15GDALRasterBandH9OGRLayerHiPPc16GDALProgressFuncPv

Fair point. I did test it and it does seem to be working on the small dataset I provided. Thanks for your help!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants