Python GDAL – Finding GeoTIFF Filename from VRT File and Coordinates

coordinatesgdalgeotiff-tiffpython

I'm new to GDAL and to GIS related Python packages in general with background more on data science libraries (such as Scikit-learn). The GIS related terminology is also a bit new to me which makes it a bit challenging to read GDAL-related documentations. Because I did not manage to find a solution with some moments of Googling or reading the GDAL manual, so I ask my question here.

My problem is very simple:

Given a virtual raster .VRT-file and some coordinate (E,N), how does one get the filename of the GeoTiff-source, which containts the given coordinate? Does there exists a function for this in GDAL?

As I've understood, .vrt-file is basically a virtual mosaic representation (stitched into one big GeoTIFF naming) of a multiple smaller GeoTiff-files (called tiles). If in Python my .vrt-file is denoted as vrt_file, then is there some function in GDAL, which allows me to give it a coordinate (E,N) and the vrt_file as input, and the function then returns the filename and/or path to the GeoTIFF, which "contains this coordinate inside its geographical coordinate borders"?

Best Answer

I solved this myself by making the following trivial script function, because I did not find a similar (and more elegant and general) from GDAL.

You give as input a coordinate of interest and a list of GeoTiff-paths which constitute the .vrt-file.

def get_geotiff_filepaths_which_contain_xy_coords(x_coord, y_coord, list_of_vrt_geotiff_filepaths):
   list_of_geotiffs_that_contain_coord_xy = []
   for input_tiff_path in list_of_vrt_geotiff_filepaths:
       tiff_data = gdal.Open(input_tiff_path, gdal.GA_ReadOnly) 
       geotransform = tiff_data.GetGeoTransform()
       # Get GeoTiff borders
       minx = geotransform[0] 
       maxy = geotransform[3] 
       pix_x = geotransform[1]
       pix_y = geotransform[5]
       x_ext = tiff_data.RasterXSize
       y_ext = tiff_data.RasterYSize 
       maxx = minx + pix_x * x_ext
       miny = maxy + pix_y * y_ext
       # Check if coordinate is inside geotiff borders
       if ((minx <= x_coord) and (x_coord <= maxx) and (miny <= y_coord) and (y_coord <= maxy)):
          list_of_geotiffs_that_contain_coord_xy.append(input_tiff_path)
   return list_of_geotiffs_that_contain_coord_xy
Related Question