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 lat_rel = site_lat - info #Get pixel size (resolution) res_we = abs(info) res_ns = abs(info) #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 + lon_rel, info, info, info + lat_rel, info, info) ds.GetRasterBand(1).WriteArray(subset) ds = None