mapwingis中的ClipGridWithPolygon方法祥读

1、祥读该方法的原因 （该方法有误）

最近在项目中用到了MapwinGIS中的

ClipGridWithPolygon (string inputGridfile, Shape poly, string resultGridfile, bool keepExtents)

和

ClipGridWithPolygon2 (Grid inputGrid, Shape poly, string resultGridfile, bool keepExtents)

 

但在我使用的时候，poly 在Grid内时没有任何问题，但 poly 的一部分在Grid外时，例如下面这种情况

就会报：

DEBUG: Image: Index Out of Bounds      错误。

查看元数据，发现只能看c++代码了。

看c中代码，你会发现 ClipGridWithPolygon (string inputGridfile, Shape poly, string resultGridfile, bool keepExtents) 方法最终还是用ClipGridWithPolygon2 (Grid inputGrid, Shape poly, string resultGridfile, bool keepExtents)了。

// **************************************************************
//		ClipGridWithPolygon()
// **************************************************************
STDMETHODIMP CUtils::ClipGridWithPolygon(BSTR inputGridfile, IShape* poly, BSTR resultGridfile, VARIANT_BOOL keepExtents, VARIANT_BOOL* retVal)
{
    AFX_MANAGE_STATE(AfxGetStaticModuleState());

    *retVal = VARIANT_FALSE;

    if (!Utility::FileExistsUnicode(inputGridfile))
    {
        ErrorMessage(tkINVALID_FILENAME);
        return S_OK;
    }

    IGrid* grid = NULL;
    CoCreateInstance(CLSID_Grid, NULL, CLSCTX_INPROC_SERVER, IID_IGrid, (void**)&grid);

    if (grid)
    {
        bool inRam = false;
        VARIANT_BOOL vb;
        grid->Open(inputGridfile, GridDataType::UnknownDataType, inRam, GridFileType::UseExtension, NULL, &vb);
        if (vb) {
            this->ClipGridWithPolygon2(grid, poly, resultGridfile, keepExtents, retVal);
            grid->Close(&vb);
        }
        grid->Release();
    }
    return S_OK;
}

 

所以继续查看ClipGridWithPolygon2 (Grid inputGrid, Shape poly, string resultGridfile, bool keepExtents)的具体代码

// ********************************************************
//		ClipGridWithPolygon2()
// ********************************************************
STDMETHODIMP CUtils::ClipGridWithPolygon2(IGrid* grid, IShape* poly, BSTR resultGridfile, VARIANT_BOOL keepExtents, VARIANT_BOOL* retVal)
{
    AFX_MANAGE_STATE(AfxGetStaticModuleState());

    *retVal = VARIANT_FALSE;

    if (!poly || !grid) {
        ErrorMessage(tkUNEXPECTED_NULL_PARAMETER);
        return S_OK;
    }

    if (Utility::FileExistsUnicode(resultGridfile))
    {
        ErrorMessage(tkFILE_EXISTS);
        return S_OK;
    }

    VARIANT_BOOL vb;

    // find cell bounds
    IExtents* ext1 = NULL;
    poly->get_Extents(&ext1);

    IExtents* ext2 = NULL;
    ((CGrid*)grid)->get_Extents(&ext2);

    IExtents* bounds = NULL;
    ((CExtents*)ext1)->GetIntersection(ext2, &bounds);

    ext1->Release();
    ext2->Release();

    if (bounds) {
        double xMin, yMin, zMin, xMax, yMax, zMax;
        bounds->GetBounds(&xMin, &yMin, &zMin, &xMax, &yMax, &zMax);
        bounds->Release();

        long firstCol, firstRow, lastCol, lastRow;
        grid->ProjToCell(xMin, yMin, &firstCol, &firstRow);
        grid->ProjToCell(xMax, yMax, &lastCol, &lastRow);

        IGrid* newGrid = NULL;

        if (!keepExtents) {
            newGrid = ((CGrid*)grid)->Clip(resultGridfile, firstCol, lastCol, firstRow, lastRow);
        }
        else {
            newGrid = ((CGrid*)grid)->Clone(resultGridfile);
        }

        // copy data
        if (newGrid)
        {
            double dx, dy, xll, yll;
            IGridHeader* header = NULL;
            newGrid->get_Header(&header);
            header->get_dX(&dx);
            header->get_dY(&dy);
            header->get_XllCenter(&xll);
            header->get_YllCenter(&yll);

            CComVariant var;
            header->get_NodataValue(&var);
            double noData;
            dVal(var, noData);

            long minRow = MIN(firstRow, lastRow);
            long maxRow = MAX(firstRow, lastRow);

            int cmnCount = lastCol - firstCol + 1;
            int rowCount = maxRow - minRow + 1;

            if (keepExtents)
            {
                long rowCount;
                header->get_NumberRows(&rowCount);

                xll += dx * firstCol;
                yll += dy * (rowCount - maxRow - 1);
            }
            header->Release();

            if (cmnCount > 0 && rowCount > 0)
            {
                double* vals = new double[cmnCount];
                long row = 0;

                CPointInPolygon pip;
                if (pip.SetPolygon(poly))
                {
                    for (long i = minRow; i <= maxRow; i++) {
                        grid->GetFloatWindow2(i, i, firstCol, lastCol, vals, &vb);

                        if (vb) {

                            //double y = yll + dy * (rowCount - row - 1.5);	// (rowCount - 1) - row;  -0.5 - center of cell
                            double y = yll + dy * (rowCount - row - 1);   // a fix suggested here: http://bugs.mapwindow.org/view.php?id=2349

                            pip.PrepareScanLine(y);

                            // set values outside polygon to nodata
                            for (long j = 0; j < cmnCount; j++)
                            {
                                // double x = xll + dx * (j + 0.5);
                                double x = xll + dx * j;		// a fix suggested here: http://bugs.mapwindow.org/view.php?id=2349

                                if (!pip.ScanPoint(x))
                                {
                                    vals[j] = noData;
                                }
                            }

                            if (keepExtents) {
                                newGrid->PutFloatWindow2(i, i, firstCol, lastCol, vals, &vb);
                                row++;
                            }
                            else {
                                newGrid->PutRow2(row++, vals, &vb);
                            }
                        }
                    }
                }
                delete[] vals;
            }

            newGrid->Save(resultGridfile, GridFileType::UseExtension, NULL, &vb);
            newGrid->Close(&vb);
            newGrid->Release();
            *retVal = VARIANT_TRUE;
        }
    }
    return S_OK;
}

2、处理方法

阅读代码你会发现源代码没什么高深的东西，要用多边形去裁切一个栅格数据，首先用多边形和栅格的外包框，去做一个相交运算。如果不相交也就没有处理的必要了。如果相交，保留原来范围只需将原grid拷贝一下，不保留原范围就是只要poly的范围，那么需要裁切grid，然后进行数据拷贝，在数据拷贝过程中，将范围外的数据设置为nodatavalue。

出现上述问题主要是poly的一部分在栅格外，而 grid->GetFloatWindow2(i, i, firstCol, lastCol, vals, &vb);却把高宽一样的数剧排除了，
 

但按逻辑修改这里后会发现RasterIO 方法又有问题，这个方法又是GDAL中的，又要去编译GDAL代码了。

 

如果要求不严可以加个判断 lastCol 和原grid的lastcol 

 

****开源虽挺好，但遇到的问题的几率也大很多，耐心阅读，慢慢修改，*****