gdalgdal-translategdalwarp
I'm new to GDAL and its utilities.
My experience with both these programs (Warp and Translate) doesn't go much further than doing some image clipping, whether by coordinate's bounding box (that I do with both programs), "size"'s bounding box (that I only can do with translate) or use a polygon cut-line (that I only can do with gdalwarp).
Besides these clipping operations, what other operations can I do? And from those which are exclusive from each program.
There are now Python modules easier to use for that, as rasterio
Rasterio employs GDAL to read and writes files using GeoTIFF and many other formats. Its API uses familiar Python and SciPy interfaces and idioms like context managers, iterators, and ndarrays.
Therefore from Masking raster with a polygon feature in Rasterio Cookbook
import rasterio
from rasterio.mask import mask
# the polygon GeoJSON geometry
geoms = [{'type': 'Polygon', 'coordinates': [[(250204.0, 141868.0), (250942.0, 141868.0), (250942.0, 141208.0), (250204.0, 141208.0), (250204.0, 141868.0)]]}]
# load the raster, mask it by the polygon and crop it
with rasterio.open("test.tif") as src:
out_image, out_transform = mask(src, geoms, crop=True)
out_meta = src.meta.copy()
# save the resulting raster
out_meta.update({"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform})
with rasterio.open("masked.tif", "w", **out_meta) as dest:
dest.write(out_image)
Result
1) You can simply use the json or the geojson modules to import the geometry
with open('poly.json') as data_file:
geoms= json.load(data_file)
print geoms
{u'type': u'Polygon', u'coordinates': [[[250204.0, 141868.0], [250942.0, 141868.0], [250942.0, 141208.0], [250204.0, 141208.0], [250204.0, 141868.0]]]}
2) you can use the ogr module with a shapefile
from osgeo import ogr
reader = ogr.Open('poly.json')
layer = reader.GetLayer()
feature = layer.GetFeature(0)
geoms =json.loads(feature.ExportToJson())['geometry']
print geoms
{u'type': u'Polygon', u'coordinates': [[[250204.0, 141868.0], [250942.0, 141868.0], [250942.0, 141208.0], [250204.0, 141208.0], [250204.0, 141868.0]]]}
3) you can also use the Fiona module
It is a Python wrapper for vector data access functions from the OGR library
import fiona
with fiona.open("poly.shp") as shapefile:
geoms = [feature["geometry"] for feature in shapefile]
print geoms
[{'type': 'Polygon', 'coordinates': [[(250204.0, 141868.0), (250942.0, 141868.0), (250942.0, 141208.0), (250204.0, 141208.0), (250204.0, 141868.0)]]}]
You can use a filter as in the script of Luke in How to set a spatial filter with Python/GDAL?. Instead of cutting, you filter the input.
from osgeo import gdal
xmin,ymin,xmax,ymax = [250204.0, 141208.0, 250942.0, 141868.0]
def map_to_pixel(mx,my,gt):
#Convert from map to pixel coordinates.
#Only works for geotransforms with no rotation.
px = int((mx - gt[0]) / gt[1]) #x pixel
py = int((my - gt[3]) / gt[5]) #y pixel
return px, py
def extent_to_offset(xmin,ymin,xmax,ymax,gt):
pxmin,pymin = map_to_pixel(xmin,ymin,gt)
pxmax,pymax = map_to_pixel(xmax,ymax,gt)
return pxmin,pymin,abs(pxmax-pxmin),abs(pymax-pymin)
ds=gdal.Open('test.tif')
gt = ds.GetGeoTransform()
bands = ds.RasterCount
band_list = []
#
# Read in bands and store all the data in bandList
#
for i in range(bands):
band = ds.GetRasterBand(i+1) # 1-based index
# apply filter to only read the data in the bounding box (xmin, ...)
data = band.ReadAsArray(xoff, yoff, xsize, ysize)
band_list.append(data)
driver = gdal.GetDriverByName("GTiff")
dst_dst = driver.Create('result.tif', xsize, ysize, 4, gdal.GDT_Byte)
for j in range(bands):
data = band_list[j]
dst_dst.GetRasterBand(j+1).WriteArray(data)
....
dst_dst = None
You only need to add the crs
You're getting a speedup using NUM_THREADS, but only at the compression stage. gdal_translate cannot used multithreading for any function apart from compression.
Probably the GDAL_CACHEMAX command is helping you out more than the NUM_THREADS option.
Best Answer
In some ways
gdalwarpis a generalization and improvement ongdal_translate.gdalwarpcan reproject rasters, this covers much of whatgdal_translatecan do (resize, crop, resample, convert format etc). It uses a different kind of settings, you specify the 'target grid properties', whereas with translate you specify 'crop/outsize'.It is also more performant and better parallelizable see arguments '-multi' and '-wo NUM_THREADS'.
https://gdal.org/programs/gdalwarp.html
'gdal_translate' is more suitable if you need to assign the extent or crs to a raster that has incorrect metadata.