0

The Situation

I use the GDAL library to warp geotiffs (subsequently referred to as "rasters") of different resolutions (30m and 500m) to a 100m (subsequently referred to as "default") resolution. The function below (WarpRasterToDefaultGridcellFormat) is supposed to be applicable to up- and downsampling tasks and should employ different aggregation algorithms:

def WarpRasterToDefaultGridcellFormat(in_rds_FilePath, defaultRaster_FilePath, polygonUniqueID, AggregationType,
                                      FileName, timer_dict=None):
    ################
    ### Link for inspiration : http://jgomezdans.github.io/gdal_notes/reprojection.html
    ### OBJECTIVE: REPROJECT A RASTER FILE SO THAT IT ALIGNS (I.E. SHARES THE SAME META-CARACTERISTICS (CRS, GEOTRANSFORM, RESOLUTION, ...)) AS A DEFAULT RASTER
    ### INPUTS: in_rds_FilePath -> PATH TO RASTER FILE TO BE REPROJECTED/up- or downsampled. TYPE: STRING
    ###         defaultRaster_FilePath -> PATH TO DEFAULT RASTER TO USE AS "TEMPLATE". TYPE: STRING
    ###         polygonUniqueID -> TILE UNIQUE ID (TILE_ID). TYPE: STRING
    ###         AggregationType -> DEFINE THE AGGREGATION TYPE. POSSIBILITIES HERE ARE SUM, AVERAGE, NEARESTNEIGHBOUR
    ###         FileName -> NAME TO BE GIVEN TO THE OUTPUT DATASET. TYPE: STRING
    ### OUTPUTS: FILEPATH TO THE CREATED RASTER DATASET.
    ################
    tic = time.perf_counter()
    in_rds = gdal.Open(in_rds_FilePath, gdalconst.GA_ReadOnly)  # GA_ReadOnly for it not to be overwritten
    in_rds_proj = osr.SpatialReference(wkt=in_rds.GetProjection())
    in_rds_geoTransform = in_rds.GetGeoTransform()

    defaultRaster_rds = gdal.Open(defaultRaster_FilePath,
                                  gdalconst.GA_ReadOnly)  # from function DefineProxyRasterDatasetBasedOnShapefileExtent
    defaultRaster_rds_proj = osr.SpatialReference(wkt=defaultRaster_rds.GetProjection())
    def_rds_geoT = defaultRaster_rds.GetGeoTransform()

    tx = osr.CoordinateTransformation(defaultRaster_rds_proj, in_rds_proj)  #

    x_size = defaultRaster_rds.RasterXSize  # number of pixels from east to west
    y_size = defaultRaster_rds.RasterYSize  #  number of pixels from north to south

    (ulx, uly, ulz) = tx.TransformPoint(def_rds_geoT[0], def_rds_geoT[3])
    # (lrx, lry, lrz) = tx.TransformPoint(def_rds_geoT[0] + def_rds_geoT[1]*x_size, def_rds_geoT[3] + def_rds_geoT[5]*y_size)   # ulx = upper left corner x, ... Function to find bottom right corner of. NOT USED?

    TileOutputFolderPath = GetTileOutputFolderPath(polygonUniqueID)
    out_rds_FilePath = os.path.join(TileOutputFolderPath, polygonUniqueID + '_' + FileName + '.tif')

    # CREATE THE OUTPUT RASTERFILE
    driver = gdal.GetDriverByName('GTiff')
    out_rds = driver.Create(out_rds_FilePath, defaultRaster_rds.RasterXSize, defaultRaster_rds.RasterYSize,
                            in_rds.RasterCount, in_rds.GetRasterBand(1).DataType, options=['COMPRESS=LZW'])

    # SET THE GEOTRANSFORM AND THE CRS OF THE OUTPUT RASTER FILE (i.e. the resampled one)
    # out_rds_geoT=(ulx, def_rds_geoT[1], def_rds_geoT[2], uly, def_rds_geoT[4], def_rds_geoT[5])
    # out_rds.SetGeoTransform(out_rds_geoT)
    out_rds.SetGeoTransform(def_rds_geoT)
    out_rds.SetProjection(defaultRaster_rds_proj.ExportToWkt())

    # REPROJECT THE INPUT RASTER USING THE APPROPRIATE AGGREGATION ALGORITHM
    if AggregationType == 'sum':
        gdal.Warp(out_rds, in_rds, resampleAlg=gdal.GRA_Sum, dstNodata=np.nan)

    elif AggregationType == 'average':
        gdal.Warp(out_rds, in_rds, resampleAlg=gdal.GRA_Average, dstNodata=np.nan)

    elif AggregationType == 'NearestNeighbour':
        gdal.Warp(out_rds, in_rds, resampleAlg=gdal.GRA_NearestNeighbour, dstNodata=np.nan)

    elif AggregationType == 'mode':
        gdal.Warp(out_rds, in_rds, resampleAlg=gdal.GRA_Mode, dstNodata=np.nan)

    # DEFINE THE APPROPRIATE RASTER DESCRIPTION
    for band_nbr in range(in_rds.RasterCount):
        band_nbr = band_nbr + 1

        RasterBand_in = in_rds.GetRasterBand(band_nbr)
        RasterBand_in_descr = RasterBand_in.GetDescription()

        RasterBand_out = out_rds.GetRasterBand(band_nbr)
        RasterBand_out.SetDescription(str(RasterBand_in_descr) + '_Reprojected')

    in_rds = defaultRaster_rds = out_rds = None

    toc = time.perf_counter()
    runtime = round((toc - tic) / 60, 2)
    print('Creation of ' + FileName + ' with the ' + AggregationType + ' algorithm took ' + str(runtime) + ' minutes.')

    if timer_dict is not None:
        timer_dict[FileName] = runtime

    return out_rds_FilePath

The problem

This approach works well for the 30m -> 100m downsampling. However, when upsampling from 500m -> 100m, the displayed resolution is 500m, but the metadata indicates that it is 100m. This is a problem regardless of the employed aggregation algorithm (e.g., NearestNeighbour, Mode, ...)

What I tried

  • Import the wrapped rasters to QGIS and check the metadata: The displayed resolution in the main window is different between the 500 -> 100m upsampled raster and the 100x100m default raster (Figure 1), although the metadata indicates that the two rasters (Figure 2 and Figure 3) should be equal.

Figure 1: Visual of the 100x100m raster, which I want to sample to (greens) and the raster that should have been upsampled to 100x100m (grayscale)

Figure 1: Visual of the 100x100m raster, which I want to sample to (greens) and the raster that should have been upsampled to 100x100m (grayscale)

Figure 2: Metadata of 100x100m default raster as a reference

Metadata of 100x100m default raster as a reference

Figure 3: Metadata of the upsampled 500 -> 100m raster, which should be equal to the metadata of the 100x100m default raster (except the statistics)

Metadata of the upsampled 500 -> 100m raster, which should be equal to the metadata of  the 100x100m default raster (except the statistics)

  • Plot the generated rasters 500 -> 100m upsampled and the 100x100m default raster, and display the metadata. Same issue, that although the metadata shows a 100x100m resolution, the 500->100m upsampled raster appears more coarse.

  • Generating the files again in Python with QGIS closed, since sometimes changes in QGIS are reflected in the files.

Ideas what the problem could be

I am pretty lost here, since the aggregation works well from 30 -> 100m and from 100m -> 10000m. The issue seems to happen only when upsampling.

MWE

Edit: Since the question is resolved, I removed the MWE from the repository.

3
  • 3
    Are you sure that the upsampled version 500->100 m just does not show equal pixel values for all the 25 pixels in the new version? Upsampling does not create new information.
    – user30184
    Jan 30 at 15:12
  • @user30184 that is it. I did not think properly that the selection of aggregation algorithms in my case does not play a role in upsampling, and therefore all the generated rasters look the same. Thanks! If you have time you could post a short answer, otherwise I can do it later with proper attribution for ideation so we can resolve the question.
    – Thierry
    Jan 30 at 16:48
  • 1
    Feel free to write an answer yourself.
    – user30184
    Jan 30 at 19:34

1 Answer 1

1

As @user30184 pointed out in the comments, upsampling does not create new information. The resolution appears more coarse since 5 adjacent pixels all have the same values.

This is by default (at least in my case) not shown through "gridlines" that would separate the single pixels, neither in QGIS nor in Python.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.