My problem is that I have multi-band images (8 bands) in GEOTIFF format and I tried to use gdal_calc.py to calculate a conditional expression between those bands. However, I could not make it because of wrong syntax. The source page does not give any example about this problem.
[GIS] How to use gdal_calc.py for multi-band images
gdalmulti-bandnumpyraster
Related Solutions
If I understand correctly you can accomplish this manually by using a tool such as Composite Bands: http://help.arcgis.com/EN/ARCGISDESKTOP/10.0/HELP/index.html#//00170000009p000000
Go into ArcCatalog and expand the multi band raster dataset (clicking the "plus" icon) and drag a single band into the input dialog for Composite Bands. Make sure the one band is the only input. Once run you will get a new single band raster dataset.
Then just repeat with the other bands.
Once you have all of the bands separated use raster to raster conversion tools to convert them to whatever format the other software you are using can consume (probably GeoTiff or ASCII). What software is it? Perhaps we could assist further in determining format.
You may find these discussion helpful:
I've tried out a bit and found that:
- RasterBricks, as mentioned by RobertH's answer, do work and are more user-friendly and easy to use;
- Rgdal methods like
readGDALalso work, but with more parameters it's a little bit less user-friendly;
So which option should one use?
According to my tests (on my 420GB GeoTiff with dimensions of 18660x21592 and 374 bands) Rgdal is faster. Maybe due to less overhead of a higher level library such as the Raster-package.
Here are my results, using system.time and replicate:
With brick:
> modis_ndvi_ts_brick <- brick("../data/pa_br_mod13q1_ndvi_250_2000_2016.tif")
> system.time(replicate(100, modis_ndvi_ts_brick[7000,7000]))
user system elapsed
94.024 5.468 99.562
While with rgdal:
> system.time(replicate(100, readGDAL("../data/pa_br_mod13q1_ndvi_250_2000_2016.tif", offset=c(7000,7000), region.dim=c(1,1))))
#some printed output from readGDAL here
user system elapsed
88.752 5.400 94.213
So as you can see Rgdal is slightly faster. For those whose this does not make a difference, I recommend using RasterBrick, it's simpler. But for those whose, like me, are struggling to create a high performance code in which every millisecond matters: Use Rgdal
Note: I'm sorry for the not exactly reproducible data, but one could try it out with other data and post here if their results differ somehow
Best Answer
To calculate a grey-scale from the same input file using different bands u can open the file multiple times and define the band which you want to use with
--A_band=n.See my example for calculating the NDVI from a satellite image with red at band 1 and near-infrared at band 4.
gdal_calc.py -A input.tif --A_band=1 -B input.tif --B_band=4 --outfile=ndvi.tif --calc="((B-A)/(B+A))"You can see that I used
input.tifforAandB. See documentation of gdal_calc.py: http://www.gdal.org/gdal_calc.html