New answers tagged

1

As stated in my comment, you should not have to cut your geotiff to bypass the 4GB limitation. There is support for Bigtiff (mentioned in the Geotiff GDAL driver documentation https://gdal.org/drivers/raster/gtiff.html) It mentioned BigTIFF is a TIFF variant which can contain more than 4GiB of data (size of classic TIFF is limited by that value). This ...


2

According to Creating geopackage from multiple raster files?, it seems you need to set RASTER_TABLE=orto_copy for the -co option instead of LAYER_NAME=orto_copy Also mentioned in https://gdal.org/drivers/raster/gpkg.html#creation-options


0

Sorry to disagree with the first answer, but using the base gdal_grid utility you can most certainly set the output resolution. You need to specify the parameter -tr ("target resolution"). See the gdal_grid man page. You would run something like: gdal_grid -tr 0.02 -a invdist:power=2 -zfield January_di <input_data> <output_interp.tif>


0

At least according to the documentation (https://docs.qgis.org/2.8/en/docs/user_manual/processing_algs/gdalogr/gdal_analysis/gridinvdist.html) there is no way to control the output raster size. EDIT: As Micha stated out in his answer, this is possible if you run the GDAL algorithm from command prompt and use -tr flag. Sorry for mislead! :) However, you can ...


3

Use the TILED=YES creation option of the GTiff driver, e.g: gdal_translate -co TILED=YES -co BLOCKXSIZE=512 -co BLOCKYSIZE=512 src.tif dst.tif


0

Handling nodata per band is always complicated whether nodata values or mask bands are used. This recent GDAL pull request seems to deal with somehow similar use case https://github.com/OSGeo/gdal/pull/3932. Be prepared to have continuous trouble with per band nodata. Test data can be generated with these commands gdal_create -outsize 100 100 -bands 3 -burn ...


0

I notice this is a rather old thread but I nevertheless have the same question I'm looking for an answer to. Would it be possible to share that script again? It seems that the dropbox link has stoped to work at this moment. Thanks, Stefan


0

Did you find a solution to this, apart from the 35000 usd maptiler pro? I am facing the same challenge with maptiles for a leaflet map where the borders between parts of the map takes so long to fix today. And would like to have an automated solution before creating bigger maps.


0

A way around. Mask raster layer > Processing > Raster values to points (Saga) Creates a shapefile with a mosaic of polygons (usually squares), one for each pixel of the mask layer. Edit mode >> select all polygons in the shapefile >> Advanced digitizing toolbar >> Merge selected features A single polygon is created that mimics exactly ...


0

Not a full solution, so you may want to iterate on it. You can use GDAL/OGR spatial SQL to do this. On the command line, you can do this easily with ogr2ogr, and dump the results to an e.g., GeoJSON file. Basically, what you want is the following code, where you need to ogr2ogr -f "GeoJSON" -dialect "SQLITE" \ -sql "SELECT f....


0

Use gdal.Warp, and pass the geometry as part of a SQL query. If you don't need to change the polygon in your Python code, just dump into a GeoJSON and read it using the cutlineDSName option for gdal.Warp. fname = '/mnt/gis_matrix/geo_tiff_1_5000_025/73634_913786_M-34-1-B-b-3-1.tif' poly = 'POLYGON (( 711263.6400000000139698 621386.0999999999767169, 711268....


0

No need to go through QGIS selection. You can go the "SQL way". There is not dependencies to QGIS: you can run it in the PyQGIS console but throught also Python script outside QGIS if GDAL Python installed from osgeo import gdal, ogr output_directory = '/tmp/' input_raster = 'OB_50M/OB_50M.tif' input_vector = 'ne_10m_admin_0_countries.shp' ...


2

You could use a "heredoc" to provide the input: gdaltransform myImage.tiff << EOF 0 0 EOF


2

You can use gdallocationinfo https://gdal.org/programs/gdallocationinfo.html but the form of the report is probably not exactly as you would like it to be. gdallocationinfo test.tif 0 0 Report: Location: (0P,0L) Band 1: Value: 32 Band 2: Value: 30 Band 3: Value: 31


0

IMO, gdal_rasterize really needs a -t_srs option to reproject vectors. I guess in traditional survey workflows vector data was often delivered in a projection suitable for direct viewing (read: any conformal projection), but these days, at least with more common "consumer-grade" datasets, it's going to be in a geographic CS about 95% of the time. ...


1

This should be enough for changing the geotransformation (captured from GDAL autotests https://github.com/OSGeo/gdal/blob/master/autotest/gcore/tiff_write.py): ds = gdal.Open('tmp/tiff57.tif', gdal.GA_Update) ds.SetGeoTransform([100, 1, 3, 200, 3, 1]) ds = None


0

Try prefixing your command with python. The default "open with" behavior for .py files is probably to display them in your editor.


0

Below is how I ended up implementing this in the Swift language. (Although my original question requested an answer for the GDAL C API, my app actually uses my own Swift language wrapper around the GDAL C functions. This abstracts the C functionality away, and makes the Swift programming much simpler. I need to add a few bits and pieces to my GDAL Swift ...


0

While I see this is an old post, I recently had the exact same issue and solved it like this: <VRTRasterBand dataType="Float32" band="1" subClass="VRTDerivedRasterBand"> <PixelFunctionType>average</PixelFunctionType> <PixelFunctionLanguage>Python</PixelFunctionLanguage> <PixelFunctionCode&...


1

It's easy using a library that's designed to handle georeferenced raster data like GDAL or rasterio (based on GDAL). Here is an example based on the docs for reprojecting with rasterio: import numpy as np import rasterio as rio from rasterio.warp import calculate_default_transform, reproject, Resampling from sklearn.impute import SimpleImputer pathhr = 'C:\\...


2

It's very easy using a library that's designed to handle georeferenced raster data like GDAL or rasterio (based on GDAL). import rasterio as rio pathhr = 'C:\\Users\\dataset\\S30W051.tif' pathout = 'C:\\Users\\dataset\\S30W051_TEST.tif' with rio.open(pathhr) as src: profile = src.profile data = src.read() data[data < 0] = 0 with rio.open(...


0

If the output is used thru web servers the best format and compression options are to make a cloud optimized geotiff COG and use oprions COMPRESS=DEFLATE PREDICTOR=YES (or 2 for integers and 3 for floats). BIGTIFF=IF_SAFER if the 4 GB file size could be surpassed.


4

You most probably do not want to compress DEM data with JPEG, that would be lossy and introduce weird steps in the data. Instead I recommend the DEFLATE compression. To improve the size savings you can also use a predictor for the compressor. See https://gdal.org/drivers/raster/gtiff.html for details and more options. for %i in (*.tif) do gdal_translate -of ...


0

You really shouldn't use the lossy JPEG algorithm for compression of DEM data. Get more info here: https://kokoalberti.com/articles/geotiff-compression-optimization-guide/ You could increase processing speed with GDAL creation options (i.e. NUM_THREADS, GDAL_CACHEMAX)


0

Apparently, this short shell script does what I want. The script is inspired from the answer that appeared in this post. #!/bin/sh for f in *.tif; do #gdal_translate -co COMPRESS=JPEG -co TILED=YES "$f" "${f%.*}.jpg" gdal_translate "$f" "jpg/${f%.*}.jpg" done


1

It looks like you're attempting to save as JPEG files. In which case GDAL would use the JPEG raster driver. The JPEG raster driver does not have a COMPRESS= creation option. See the driver's documentation at https://gdal.org/drivers/raster/jpeg.html#raster-jpeg and scroll down to see the list of creation options. In fact, setting a compression option to ...


5

You can first convert all tif input as a single vrt, then run the gdal.Warp from pathlib import Path from osgeo import gdal vrt_name= 'input.vrt' tifs = [str(p) for p in Path('.').glob('*.tiff')] my_vrt = gdal.BuildVRT(vrt_name, tifs) my_vrt = None ds_input = gdal.OpenEx(vrt_name, gdal.OF_RASTER) ds = gdal.Warp('output.tif', ds_input, options="-...


1

Turns out I needed to use gdalbuildvrt -resolution highest


2

You can greatly simplify the operation using gdal.VectorTranslate from osgeo import gdal # Open the raster tif and extract crs t = "path/to/raster.tif" tif = gdal.Open(t) crs1 = tif.GetProjectionRef() tif = None # Input shp path s = "path/to/shapefile.shp" # Set name of the new shapefile sUTM = "shapefile_UTM.shp" if os.path....


1

I'll add this as an answer (was originally meant as a comment). Perhaps it helps someone since there is and has been a lot of confusion about this. It expands the answers above. There is another helpful question with answers. In short, to print the wkt comment and avoid the warnings with PROJ6 CRS objects, use sp::wkt instead of sp::proj4string: sp::wkt(x) [...


0

The problem has been resolved. The point was that when creating the original Geotiff, I designated nodata = 0, which was an error, since I work with NDVI. Changing the nodata value in the source file to any other solves the problem.


2

The cause: It seems to be due to a PostgreSQL change that occurred in version 12.0. Namely, the extra_float_digits parameter that is now being used for pg_dump and pg_dumpall functions (see here: https://www.postgresql.org/docs/release/12.0/ and here: https://www.postgresql.org/docs/12/runtime-config-client.html#GUC-EXTRA-FLOAT-DIGITS). The solution: Set the ...


0

You need to transform your data in .tiff with gdal_translate, make the resize with gdalwarp, and then back to .nc with gdal_translate again. Tips for increased speed: https://trac.osgeo.org/gdal/wiki/UserDocs/GdalWarp


0

The idiomatic algorithm for extraction of pixels closest to a line is the Bresenham algorithm. It may be present in underlying guts of libraries such as GDAL or others. You will most probably find lots of implementations in various languages with some web search. Quoting wikipedia: Bresenham's line algorithm is a line drawing algorithm that determines the ...


1

If the aim is just to divide the image into two values, one for data and another for nodata and exact values are not so important, a short gdal_calc https://gdal.org/programs/gdal_calc.html command can do it. gdal_calc -A src.tif --outfile=out.tif --type=Byte --calc="255*(A>-9999)" The command creates a single band, 8-bit GeoTIFF because there ...


2

I am not sure if this is a GDAL issue. Perhaps it is not an issue at all but rounding errors just belong to conversions between decimal numbers and binary values that computers are using. I stored your data into PostGIS with GDAL 3.3.0 and the geometry appears to be in the database with unaltered precision when queried this way: select ST_AsText(wkb_geometry)...


Top 50 recent answers are included