Clipping area from big raster without Resampling

Question

I'm trying to polygonize a raster. The raster is a very big raster so I try to clip it with gdalwrap and a polygon in a shp file as mask.

After cliping (and polygonize it) I see that the pixel are moved and in some cases the value is changed (following resampling_method parameter).

What I want is that the pixels are in the same place that the original one, without offset and without re sampling, im not worry about the final raster had more (not less) area that mask polygon.

The trasparent area is the mask (the area were I want to work). The first image is the poligonized result without cliping. The second one is cliped in a bigger poligon (pixel still moved and modified). The third one is cliped to the interes area with this command: gdalwarp --config GDALWARP_IGNORE_BAD_CUTLINE YES -r max -q -crop_to_cutline -cutline $3 -of GTiff raster.TIF maskFile

With a bigger mask area

If you take a close look to the square you will see that the pixels are moved, and values changed because pixels with the same value are joined.

What I want is something like second image but with the pixels origins and values equals to original raster. (Not worried if I have to do in code or command line or wherever)

Answer 1

As I get no answers I did my own code, I put it here for it someone want it.

#include "gdal/gdal_priv.h"
#include "gdal/cpl_conv.h" // for CPLMalloc()
#include "gdal.h"
#include "gdal_alg.h"
#include "ogr_srs_api.h"
#include "cpl_string.h"
#include "cpl_conv.h"
#include <ogrsf_frmts.h>
#include <stdlib.h>
#include <dirent.h>
#include <list>



void getShapeLimits(GDALDataset* shpDS,OGREnvelope *mL ) {
    OGRLayer *ly1;
    int elements;
    ly1=shpDS->GetLayer(0);
    elements=ly1->GetFeatureCount();
    ly1->GetExtent(mL);
}

void getRasterLimits(GDALDataset*originDS,OGREnvelope *mL ) {
    int cols = originDS->GetRasterXSize();
    int rows = originDS->GetRasterYSize();
    double originTransform[6];
    originDS->GetGeoTransform(originTransform);
    mL->MinX = originTransform[0];
    mL->MaxX = originTransform[0] + cols * originTransform[1] + rows * originTransform[2];

    mL->MinY = originTransform[3] + cols * originTransform[4] + rows * originTransform[5];
    mL->MaxY = originTransform[3];


}
bool checkInclusion(GDALDataset*areaDS,OGREnvelope &shapeEV)
{
    int cols = areaDS->GetRasterXSize();
    int rows = areaDS->GetRasterYSize();
    if((shapeEV.MinX <0) || (shapeEV.MaxX)>cols)
        return 0;
    if((shapeEV.MinY>rows)||(shapeEV.MaxY<0))
        return 0;
    return 1;
}

int main(int argc, char *argv[])
{
    GDALDataset *origDS,*shapeDS,*destinationDS;
    double shapeTransform[6];
    double originTransform[6];
    OGRRegisterAll();
    if(argc<4) {
        printf("Error: landsatReduceC origin.tiff interesArea.shp output.tiff");
        fflush(stdout);
        exit(-1);
    }
    origDS =  (GDALDataset*) GDALOpen(argv[1],GA_ReadOnly);
    if(!origDS) {
        printf( "Error opening %s\n", argv[1] );
        exit( -1 );
    }
    shapeDS =  (GDALDataset*) GDALOpenEx( argv[2], GDAL_OF_VECTOR, NULL, NULL, NULL );
    if(!shapeDS) {
        printf( "Error opening %s\n", argv[2] );
        exit( -1);
    }
    OGREnvelope shapeLimits,originLimits;
    getRasterLimits(origDS,&originLimits);

    getShapeLimits(shapeDS,&shapeLimits);

    origDS->GetGeoTransform(originTransform);
    shapeDS->GetGeoTransform(shapeTransform);
    OGRSpatialReference *oSourceSRS, oTargetSRS;
    OGRErr error;
    OGRCoordinateTransformation *poCT;
    char *pszWkt=(char *)origDS->GetProjectionRef();
    error = oTargetSRS.importFromWkt(&pszWkt);
    oSourceSRS=shapeDS->GetLayer(0)->GetSpatialRef();
    error = oSourceSRS->exportToWkt(&pszWkt);
    poCT = OGRCreateCoordinateTransformation(oSourceSRS,&oTargetSRS);

    double x[2];
    double y[2];
    x[0]=shapeLimits.MinX;
    y[0]=shapeLimits.MaxY;
    x[1]=shapeLimits.MaxX;
    y[1]=shapeLimits.MinY;

    poCT->Transform(2,x,y);
    OGREnvelope clipArea;
    int secMarginMts = 1;
    int xtl = clipArea.MinX = (x[0] - originTransform[0]) / originTransform[1]-secMarginMts;
    int ytl = clipArea.MaxY = (y[0] - originTransform[3]) / originTransform[5]-secMarginMts;
    int xbr = clipArea.MaxX = (x[1] - originTransform[0]) / originTransform[1]+secMarginMts;
    int ybr = clipArea.MinY = (y[1] - originTransform[3]) / originTransform[5]+secMarginMts;

    if(!checkInclusion(origDS,clipArea)) {
        printf("Out of limits\r\n");
        exit(-1);
    }

    GDALRasterBand *originalRB = origDS->GetRasterBand(1);
    GDALDataType dt=originalRB->GetRasterDataType();
    if(dt!=GDT_UInt16) {
        printf("uint16 raster expected and we get type%d (check GDALDataType)\r\n",dt);
        exit(-1);
    }
    uint16_t *clipedData;
    clipedData=(uint16_t *)malloc(sizeof(uint16_t)*(xbr-xtl)*(ybr-ytl));
    originalRB->RasterIO(GF_Read,xtl,ytl,xbr-xtl,ybr-ytl,clipedData,xbr-xtl,ybr-ytl,dt,0,0);

    if(clipedData[0]==1){
        printf("This image has no info");
        exit(-1);
    }
    //make the raster
    double new_x = originTransform[0] + xtl*originTransform[1];
    double new_y = originTransform[3] + ytl*originTransform[5];

    double clipedTransform[6];
    clipedTransform[0]=new_x;
    clipedTransform[1]=originTransform[1];
    clipedTransform[2]=originTransform[2];
    clipedTransform[3]=new_y;
    clipedTransform[4]=originTransform[4];
    clipedTransform[5]=originTransform[5];
    GDALDriver *poDriver;
    poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
    destinationDS = poDriver->Create(argv[3],xbr-xtl,ybr-ytl,1,GDT_UInt16,NULL);
    destinationDS->GetRasterBand(1)->RasterIO(GF_Write,0,0,xbr-xtl,ybr-ytl,clipedData,xbr-xtl,ybr-ytl,GDT_UInt16,0,0);

    destinationDS->SetGeoTransform(clipedTransform);
    destinationDS->SetProjection(origDS->GetProjectionRef());

    GDALClose(destinationDS);
    exit(0);


}