0

I have a TIFF file which I need to convert to dataframe and back to TIFF. This is what I did so far:

from osgeo import gdal
import pandas as pd
import numpy as np

ds = gdal.Open('t2m_20210523.tif')

ar = ds.GetRasterBand(1).ReadAsArray()
flat = ar.flatten()
gt = ds.GetGeoTransform()
res = gt[1]
xmin = gt[0]
ymax = gt[3]
xsize = ds.RasterXSize
ysize = ds.RasterYSize
xstart = xmin +res/2
ystart = ymax - res/2
ds = None
x = np.arange(xstart, xstart+xsize*res, res)
y = np.arange(ystart, ystart-ysize*res, -res)
x = np.tile(x, ysize)
y = np.repeat(y, xsize)

dfn = pd.DataFrame({"x":x[0:len(flat)], "y":y[0:len(flat)], "value":flat})
dfn['value']=dfn['value'].fillna(-999)
dfn=dfn.sort_values(['y','x'])

dfn.to_csv("df.xyz", index = False, header = None, sep = " ")
demn = gdal.Translate("out_file.tif", "df.xyz", outputSRS = "EPSG:4326")
demn = None

This is what I got: enter image description here I believe I somehow flipped coordinates, but so far I failed to find combination which works.

2
  • 1
    Why go through CSV? If your data is already gridded then try creating a new dataset then using ds.WriteArray() to populate it. Or better yet, do the same using rasterio.
    – mikewatt
    May 20, 2021 at 23:21
  • Can't you just use e.g. gdal.Translate on the original image with option outputBounds to do the half-pixel shift? You're in effect just changing the geotransform, not interpolating the data as far as I can tell.
    – Jose
    May 21, 2021 at 13:44

1 Answer 1

4

From your code it appears that you're trying to move the image half a pixel in the x and y direction. As @mikewatt suggests above this is much easier using rasterio and affine.

The affine library is a simple wrapper around Affine transformations. This is important when you're dealing with raster data as the coordinates of raster data consist of:

  1. The map projection (which describes the measurement units among other things)
  2. The GeoTransform which are the parameters for an affine transformation which translate between projected units and pixel index (ie. xy to ij)

You want to shift half a cell, and happily it's easy to apply an affine transformation to another transformation. So your first step is to create a new affine transformation.

import rasterio
from affine import Affine

src = rasterio.open('t2m_20210523.tif')


# src.transform is the affine transform for the source raster (the GeoTransform)
x_cellsize = src.transform.a
y_cellsize = src.transform.e  # will be negative

# a simple shift by half a cell
shift = Affine(1, 0, x_cellsize / 2, 0, 1, y_cellsize / 2)
new_transform = src.transform * shift

At this point it's a simple matter to write the data back out to disk:

new_meta = src.meta.copy()
new_meta['transform'] = new_transform

with rasterio.open('out_file.tif', 'w', **new_meta) as dest:
    dest.write(src.read())

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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