6

I'm sub-setting a GeoTIFF with the GDAL Python binding and numpy. Essentially, I am providing a central location (coordinates of a flux tower) from which I produce a subsetted image. I then want to write this to a new GeoTIFF. However, my current code does not correctly georeference. I'm having a hard time digesting the rather sparse and confusing GDAL documentation and so would love for some help. I have read that this might be possible with a gdal_translate or a VRT, but can not figure out any more from the docs.

So, here is how I load the GeoTIFF in question and produce the subset:

datafile = gdal.Open("%s" % filenames_matched[site][filen] )
data = datafile.ReadAsArray()
info = datafile.GetGeoTransform()       
# GDAL provides the top-left coordinate (lat, lon)
#Site_lon and site_lat are the coordinates for the fluxtower site.
# Produce a relative coordinate for the fluxtower (eg relative to gdal ul(lat,lon))
lon_rel = site_lon - info[0]
lat_rel = site_lat - info[3]

#Get pixel size (resolution)
res_we = abs(info[1])
res_ns = abs(info[5])   

#Now convert the relative position from degrees to pixels..
pixel_lon = float ( lon_rel) / res_we
pixel_lat = float ( lat_rel ) / res_ns

# Will, of course, have to round this to the nearest whole pixel...
pixel_lon = round( pixel_lon )
pixel_lat = round( pixel_lat )

# Now need to outline the actual buffer of 24 on each side...
buffer_lon_east = pixel_lon + 24.5
buffer_lon_west = pixel_lon - 24.5
buffer_lat_north = pixel_lat + 24.5
buffer_lat_south = pixel_lat - 24.5

#Produce subset array
subset = data[buffer_lon_west:buffer_lon_east,buffer_lat_south:buffer_lat_north]

And then save the new subsetted GeoTIFF:

Nx, Ny = subset.shape
driver = gdal.GetDriverByName("GTiff")

projection = datafile.GetProjection()    # Same projection as input (WGS84)
ds = driver.Create("%s" % composited_filename, Nx, Ny, 1, gdal.GDT_Float32)
ds.SetProjection(projection)

# Where I imagine the problem is:
# Sets new coordinates from derived ul lat, lon for the subset 
ds.SetGeoTransform([info[0] + lon_rel, info[1], info[2], info[3] + lat_rel, info[4], info[5])

ds.GetRasterBand(1).WriteArray(subset)
ds = None
5

You can indeed use gdal_translate either giving the source coordinates ([-srcwin xoff yoff xsize ysize]) or the georeferenced coordinates ([-projwin ulx uly lrx lry]). For instance:

import os
inDS = ... # input raster
outDS = ... # output raster
lon = ... # lon of your flux tower
lat = ... # lat of your flux tower
ulx = lon - 24.5
uly = lat + 24.5
lrx = lon + 24.5
lry = lat - 24.5
translate = 'gdal_translate -projwin %s %s %s %s %s %s' %(ulx, uly, lrx, lry, inDS, outDS)
os.system(translate)
1

You can use gdal_translate to solve this problem. However if you wanted to avoid using GDAL command line or making system calls out to it in Python, you can do it using the Python GDAL bindings directly doing something like this:

def subsection(file_name, someNewMinX, someNewMinY, someNewMaxX, someNewMaxY):
    raw_file_name = os.path.splitext(os.path.basename(file_name))[0]
    driver = gdal.GetDriverByName('GTiff')
    dataset = gdal.Open(file_name)
    band = dataset.GetRasterBand(1)
    transform = dataset.GetGeoTransform()

    cols = dataset.RasterXSize
    rows = dataset.RasterYSize
    transform = dataset.GetGeoTransform()
    minx = transform[0]
    maxx = transform[0] + cols * transform[1] + rows * transform[2]

    miny = transform[3] + cols * transform[4] + rows * transform[5]
    maxy = transform[3]

    width = maxx - minx
    height = maxy - miny

    output_path = os.path.join("data", raw_file_name)
    if not os.path.exists(output_path):
        os.makedirs(output_path)

    tiles = create_tiles(minx, miny, maxx, maxy, n)
    transform = dataset.GetGeoTransform()
    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = -transform[5]

    print xOrigin, yOrigin

    # Subsitute with your new subsection values
    newminx = someNewMinX
    newmaxx = someNewMaxX
    newminy = someNewMinY
    newmaxy = someNewMaxY

    p1 = (newminx, newmaxy)
    p2 = (newmaxx, newminy)

    i1 = int((p1[0] - xOrigin) / pixelWidth)
    j1 = int((yOrigin - p1[1])  / pixelHeight)
    i2 = int((p2[0] - xOrigin) / pixelWidth)
    j2 = int((yOrigin - p2[1]) / pixelHeight)

    print i1, j1
    print i2, j2

    new_cols = i2-i1
    new_rows = j2-j1

    data = band.ReadAsArray(i1, j1, new_cols, new_rows)

    #print data

    new_x = xOrigin + i1*pixelWidth
    new_y = yOrigin - j1*pixelHeight

    print new_x, new_y

    new_transform = (new_x, transform[1], transform[2], new_y, transform[4], transform[5])

    output_file_base = raw_file_name + "_" + "subsection" + ".tif"
    output_file = os.path.join("data", raw_file_name, output_file_base)

    dst_ds = driver.Create(output_file,
                           new_cols,
                           new_rows,
                           1,
                           gdal.GDT_Float32)

    #writting output raster
    dst_ds.GetRasterBand(1).WriteArray( data )

    #setting extension of output raster
    # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
    dst_ds.SetGeoTransform(new_transform)

    wkt = dataset.GetProjection()

    # setting spatial reference of output raster
    srs = osr.SpatialReference()
    srs.ImportFromWkt(wkt)
    dst_ds.SetProjection( srs.ExportToWkt() )

    #Close output raster dataset
    dst_ds = None

    dataset = None

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.