38
It's a well known and longstanding issue that gdalwarp doesn't deal with compression well. The solution is to gdalwarp without compression then gdal_translate with compression.
To avoid two lengthy processes, gdalwarp to VRT first, it's really quick, then gdal_translate with the -co compress=lzw option.
i.e.
$ gdalwarp -tap -tr 30 30 -t_srs "etc..." -of ...
29
Using rasterio you could do
import rasterio
file_list = ['file1.tif', 'file2.tif', 'file3.tif']
# Read metadata of first file
with rasterio.open(file_list[0]) as src0:
meta = src0.meta
# Update meta to reflect the number of layers
meta.update(count = len(file_list))
# Read each layer and write it to stack
with rasterio.open('stack.tif', 'w', **meta) ...
28
I propose two solutions: the first one using QGIS, the second one using Python (GDAL).
Solution using QGIS
In QGIS you may create a VRT mosaic.
Please follow this procedure (see the image below):
Load the raster in the Layers Panel;
Right-click on it and choose Save As...;
Check the Create VRT option;
Choose the folder where your outputs will be saved;
...
22
Your script is missing the ds.FlushCache method, that saves to disk what you have in memory at the end of the modifications. See below a corrected version of your example. Notice that I also added two lines to set projection and geotransform as input
import os
import gdal
file = "path+filename"
ds = gdal.Open(file)
band = ds.GetRasterBand(1)
arr = ...
19
So you want a geotiff instead of a tiff with a world file (.tfw). This should be the default in GDAL (http://www.gdal.org/frmt_gtiff.html) so:
gdal_translate -of GTiff input.tif ouput.tif
It will default to geotiff.
19
I would recommend to use gdalcopyproj.py, a sample file from the GDAL repository done for this purpose as mentioned directly in the script:
Duplicate the geotransform and projection metadata from
one raster dataset to another, which can be useful after
performing image manipulations with other software that
ignores or discards georeferencing ...
19
This should get you going. The raster values are read using rasterio, and pixel centre coordinates are converted to Eastings/Northings using affine, which are then converted to Latitude/Longitude using pyproj. Most arrays have the same shape as the input raster.
import rasterio
import numpy as np
from affine import Affine
from pyproj import Proj, transform
...
18
If QGIS is runnig in a 1000x1000 pixel sized window on your screen there is no need to read all 32000x32000 pixels for showing the map. GDAL tries to read data from the source image so that no data at all is read outsize the bounding box, and if image has overviews the data come from the resolution level that is best suitable for the map resolution. There is ...
18
You can do this using GDAL, it directly supports XYZ format. It doesn't matter if your coordinates are UTM, gdal_translate will output in the same coordinate system.
So to convert to GeoTIFF is as simple as:
gdal_translate test.xyz test.tif
Look at the GeoTIFF doc for output options (such as compression) and the gdal_translate doc for more usage info. In ...
17
Building on what @David mentioned you may use open source gdal library using python module to get image extent like this:
import gdal
from gdalconst import GA_ReadOnly
data = gdal.Open('C:/Temp/myimage.tif', GA_ReadOnly)
geoTransform = data.GetGeoTransform()
minx = geoTransform[0]
maxy = geoTransform[3]
maxx = minx + geoTransform[1] * data.RasterXSize
miny ...
15
If using GDAL 2.1+ it's as simple as gdal.BuildVRT then gdal.Translate:
from osgeo import gdal
outvrt = '/vsimem/stacked.vrt' #/vsimem is special in-memory virtual "directory"
outtif = '/tmp/stacked.tif'
tifs = ['a.tif', 'b.tif', 'c.tif', 'd.tif']
#or for all tifs in a dir
#import glob
#tifs = glob.glob('dir/*.tif')
outds = gdal.BuildVRT(outvrt, tifs, ...
15
Create a box as a Spatial object and crop your raster by the box.
e <- as(extent(-16, -7.25, 4, 12.75), 'SpatialPolygons')
crs(e) <- "+proj=longlat +datum=WGS84 +no_defs"
r <- crop(worldpopcount, e)
14
Short answer: I suspect there is no such standard for TIFF or GeoTIFF overviews. There are multiple implementations, methods and formats to define overviews for GeoTIFFs.
GeoTIFF is based on the TIFF format (PDF specification for Revision 6.0, from 1992). The format has support for multi-page documents or subfiles, similar to a multi-page PDF.
A GeoTIFF ...
14
I haven't found any specific commandline utility that can report if a tiff is tiled or striped. At least not directly or in a grepable form like TILED=YES.
There should be enough information in gdalinfo to make that decision, however.
I have a landsat scene, each made with gdal_translate:
landsat_tiled.tif : -co TILED=YES
landsat_notiled.tif: -co TILED=...
14
New raster tables must be added into existing GeoPackage as subdatasets as documented in the GDAL driver manual page http://www.gdal.org/drv_geopackage_raster.html. Have a look at the fifth example
Addition of a new subdataset to an existing GeoPackage, and choose a
non default name for the raster table.
gdal_translate -of GPKG new.tif existing.gpkg -...
13
gdal_translate -of GTiff C:\temp\input\a.img C:\temp\output\a.tif
and the batch option
GDAL_translate: converting ESRI GRID to Geotiff in batch
'-of GTiff' this part is probably not even required as this is the default but will do no harm.
13
You can use Rasterio to get the bounding box as follows:
import rasterio
dataset = rasterio.open('example.tif')
>>> dataset.bounds
BoundingBox(left=358485.0, bottom=4028985.0, right=590415.0, top=4265115.0)
13
Are there any drawbacks to doing this?
There are softwares that are not going to be able to read BIGTIFF. It is likely that most geospatial/GIS software that consumes raster data can now handle BIGTIFF, but others might not have been updated or make assumptions about maximum sizes. BIGTIFF support did and does require some software implementation.
If you ...
12
Quantum GIS now supports this feature, the resolution can be set and can include an optional world file containing georeferencing information.
In the print composer check the "World file on" box under Export settings.
12
The coordinates of the target extent have to be expressed in the target SRS:
-te xmin ymin xmax ymax:
set georeferenced extents of output file to be created (in target SRS).
Being...
>cs2cs +init=EPSG:4326 +to +init=EPSG:3857
5 43
556597.45 5311971.85 0.00
15 48
1669792.36 6106854.83 0.00
the command should be something like:
...
12
You have two options here:
1) Create a gdal.WarpOptions object and pass it to gdal.Warp as the options argument, just as you are doing (although you can skip the fourth line as you are creating the same object on the fifth line).
2) Create a dictionary with the arguments and pass it as keyword arguments to gdal.Warp. For example:
kwargs = {'format': '...
answered Nov 13 '19 at 14:05
Marcelo Villa-Piñeros
5,5882 gold badges15 silver badges34 bronze badges
12
Assuming your rasters are integer-type, you can create a ColorTable, specify the color for each value using the SetColorEntry() method and then apply the ColorTable to the raster using the SetRasterColorTable() method to the individual band. The SetColorEntry(pixel_val, (r, g, b)) method takes two arguments, where the first one is the pixel value and the ...
answered Jun 13 '19 at 16:44
Marcelo Villa-Piñeros
5,5882 gold badges15 silver badges34 bronze badges
11
here's a better answer, use gdalbuildvrt with either srcnodata or vrtnodata flag:
gdalbuildvrt -srcnodata "123 231 67" outfile.vrt input.tif
If the next application in line doesn't understand .vrt, translate to a new tif:
gdal_translate outfile.vrt final.tif
11
You can easily open ASCII xyz triplicate data in QGIS under "Add Raster Data" with a "ASCII Gridded XYZ (.xyz)" file type. You can also covert it to a different format under the "Raster > Conversion > Translate (Convert format)" menu. Alternately, you can do this under the "Raster > Conversion > Rasterize" menu with a "Comma Separated Value (.csv)" file type....
11
Those files are generated in ArcGIS. Are auxiliary files for georeferencing and visualization, aren't part of GTiff as .shx or .dbf are in a shapefile.
.tfw is ESRI World file
.ovr are piramid layers
.xml is schema look and histogram
.cpg is for TIFF interpretation.
.dbf is for raster attribute table (thanks to @radouxju)
You can add .tfw using GDAL ...
10
Or you can use geotifcp (http://www.remotesensing.org/geotiff/geotifcp.html).
To dump the data from a GeoTiff to a world file try listgeo (http://www.remotesensing.org/geotiff/listgeo.html).
Updated links as of 2017-05-21
http://geotiff.maptools.org/geotifcp.html
http://geotiff.maptools.org/listgeo.html
10
If you have a raster DEM already, then there is a tool that I developed in Whitebox Geospatial Analysis Tools called Remove Off-Terrain Objects, contained within the LiDAR toolbox, that works well for creating bare-earth DEMs, particularly in urban and agricultural settings. It works less well where either the terrain is steeply sloped or the forest cover is ...
10
There are two programs using the libgeotiff library that will let you export and import georeferencing info from GeoTIFFs (and other image file formats): listgeo for saving the georeferencing info into a file, and geotifcp for writing the info from that file back to the image file after it's been edited. Downloads for the libgeotiff library can be found here:...
10
Based on the advice in "GeoServer on Steroids" I would aim for a mosaic of GeoTiffs. Page 8 clearly states to choose a mosaic when:
A single file gets too big (inefficient seeks, too much metadata
to read, etc..)
Single granules can be large
Use Tiling + Overviews + Compression on granules
but suggests moving to pyramids when
...
10
If you read a TIFF with raster like you have in your loop you'll only get one layer read:
> r = raster("./GcrfPicture.tif")
> r
class : RasterLayer
band : 1 (of 4 bands)
dimensions : 720, 960, 691200 (nrow, ncol, ncell)
resolution : 1, 1 (x, y)
The 1 (of 4 bands) is telling you that the source had four bands but its only read one....
Only top voted, non community-wiki answers of a minimum length are eligible