Can be done fairly easy. First, create the dataset:
drv = gdal.GetDriverByName("GTiff")
ds = drv.Create("name.tif", width, height, 6, gdal.GDT_Float32)
You then need to set the geotransform using ds.SetGeoTransform where the argument is a six element tuple: (upper_left_x, x_resolution, x_skew, upper_left_y, y_skew, y_resolution)
Set the projection using ds.SetProjection, argument is the way string, then
Write the array for the first band:
ds.GetRasterBand(1).WriteArray(arr)
That should be enough to get you going. If the purpose is just to reproject the data, then you could look at using an in memory dataset
I dont't know, but it's easy to do with osgeo.gdal, GeoPandas and shapely box ( shapely.geometry.box(minx, miny, maxx, maxy, ccw=True) = bounding box)ΒΆ
Original raster files
import os
StartDir = "/Shared/scan_ign/68"
for dir, subdir, files in os.walk(StartDir):
for fname in files:
if fname.endswith(".tif"):
print(fname,)
68_1.tif
68_2.tif
68_3.tif
68_4.tif
68_5.tif
68_6.tif
68_7.tif
68_8.tif
With gdalindex
gdaltindex index.shp *.tif
Control
import geopandas as gpd
index = gpd.read_file('index.shp')
print index.head()
location geometry
0 68_1.tif POLYGON ((226018.4754020752 58174.93482261647,...
1 68_2.tif POLYGON ((234019.0396324735 58174.8250049785, ...
2 68_3.tif POLYGON ((242020.0156871924 58175.55617328714,...
3 68_4.tif POLYGON ((250020.4746237129 58176.34366682862,...
4 68_5.tif POLYGON ((226019.5264973617 48175.2918232924, ...
With osgeo.gdal
from osgeo import gdal
import geopandas as gpd
from shapely.geometry import box
# compute the bounding box of a gdal raster file
def bounds_raster(path):
raster = gdal.Open(path)
ulx, xres, xskew, uly, yskew, yres = raster.GetGeoTransform()
lrx = ulx + (raster.RasterXSize * xres)
lry = uly + (raster.RasterYSize * yres)
return box(lrx,lry,ulx,uly)
# creation of the index file
df = gpd.GeoDataFrame(columns=['location','geometry'])
# iterate through multiple tif files in a folder
for dir, subdir, files in os.walk(StartDir):
for fname in files:
if fname.endswith(".tif"):
df = df.append({'location':fname, 'geometry': bounds( os.path.join(dir+"/", fname))},ignore_index=True)
print df.head()
location geometry
0 68_1.tif POLYGON ((226018.4754020752 48173.94230128427,...
1 68_2.tif POLYGON ((234019.0396324735 48173.84503544428,...
2 68_3.tif POLYGON ((242020.0156871924 48174.31088528392,...
3 68_4.tif POLYGON ((250020.4746237129 48174.72574336932,...
4 68_5.tif POLYGON ((226019.5264973617 38173.0677178388, ...
# save resulting shapefile
df.to_file("tile-index.shp")
As Spacedman points out, you can also use rasterio
import rasterio
df = gpd.GeoDataFrame(columns=['location','geometry'])
for dir, subdir, files in os.walk(StartDir):
for fname in files:
if fname.endswith(".tif"):
bounds =rasterio.open(os.path.join(dir+"/", fname)).bounds
df = df.append({'location':fname, 'geometry': box(bounds[0], bounds[1], bounds[2], bounds[3])},ignore_index=True)
df.to_file("tile-index2.shp")
Best Answer
You can use an optfile (see the bottom of the page) which contains paths to all your tif files.
You can create the optfile using a dos command prompt - something like:
You can then just call gdaltindex with the optfile:
As for the different projections, if you know what projectios everything is in you can use a different optfile. Otherwise you could potentially do the following (bit sloppy I know, but I can't think of a way you could do it otherwise):