42

You need to use the -dstalpha option to gdalwarp e.g.: gdalwarp -cutline INPUT.shp -crop_to_cutline -dstalpha INPUT.tif OUTPUT.tif This will add an alpha band to the output tiff which masks out the area falling outside the cutline. A late answer, but hopefully it will help someone else with the same problem.


33

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 ...


18

Here is an example that does roughly what you ask for. The main parameters are the geotransform array that gdal uses to describe a raster location (position, pixel scale, and skew) and the epsg code of the projection. With that, the following code should properly georeference the raster and specify its projection. I did not test this much, but it seemed to ...


17

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 ...


17

I get the same results as gdalwarp from gdal.AutoCreateWarpedVRT if I set the error threshold to 0.125 to match the default (-et) in gdalwarp. Alternatively, you could set -et 0.0 in your call to gdalwarp to match the default in gdal.AutoCreateWarpedVRT. Example Create a reference to compare to: gdalwarp -t_srs EPSG:4326 byte.tif warp_ref.tif Run the ...


15

I just happened to come across this question and a potential answer when looking for something else. gdal_merge.py uses nearest neighbor resampling. If you want control over the resampling used, you should use gdalwarp instead. source: trac.osgeo.org


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

if you want to set the projection of a file with gdal tools, you should use gdal_edit.py -a_srs Both gdalwarp and gdal_translate will create a new file.


11

I'm agree with Nathan. You need to pythonize your whole script. So substitute your for loop with something like the following: import os, fnmatch def findRasters (path, filter): for root, dirs, files in os.walk(path): for file in fnmatch.filter(files, filter): yield file for raster in findRasters(INPUT_FOLDER, '*.tif'): ...


11

Reorganise your shapefile so that one shapefile contains one feature (A,B,C in your case) only Then use a loop like for i in A B C; do gdalwarp -cutline $i.shp ... $i.tif done to create each output raster. Example of script: #!/bin/sh # "shp" - folder for shapefiles # "outputraster" - folder for output rasters cd /home/user/cliprasters/ # ...


11

I don't know if it's possible to clip a raster with an other raster but you could use gdaltindex to build the shapefile with the extent of your raster. http://www.gdal.org/gdaltindex.html


11

http://www.gdal.org/gdalwarp.html says -t_srs srs_def: target spatial reference set. The coordinate systems that can be passed are anything supported by the OGRSpatialReference.SetFromUserInput() call, which includes EPSG PCS and GCSes (i.e. EPSG:4296), PROJ.4 declarations (as above), or the name of a .prj file containing well known text. ...


11

gdalwarp is doing the right thing: preserving total resolution of your image by changing the pixel-size. WGS 84 / Pseudo-Mercator projection is heavily distorted when moving away from the equator. Thus, it could be discussed if the units should be called "Pseudo-meters". One meter in reality is approximately 1/cos(lat) pseudo-meters. You can calculate the ...


10

Nice and reproducible question. Personally, I'd expect that the reason for the difference is in the implementations of the bilinear reprojection. You can obviously look into source code for the two approaches, but I'd expect that to be a vast overkill. It appears that the R implementation introduces bigger "errors" / "changes" than the raw GDAL version (...


9

Firstly, crop the image source (coords are expressed in pixels here) with: gdal_translate -srcwin 115 18 1360 2156 2104.gif 2104_cropped.tif Then, transform the known WGS84 coordinates of the upper left and lower right corners to the "WGS 84 / World Mercator" projection (EPSG:3395): cs2cs +init=epsg:4326 +to +init=epsg:3395 -35 75 -3896182.18 12890914....


8

I noticed that gdal2tiles numbers the tiles from south to north (according to the TMS specification), while Openstreetmap and others do it from north to south. For my personal use, I changed the code of gdal2tiles to get it right again. See also: http://osgeo-org.1560.x6.nabble.com/gdal2tiles-tiles-in-wrong-hemisphere-and-or-Openlayers-problem-td3742809....


8

For irregular polygons, and assuming that your geotiff raster file is a binary raster, you could use GDAL_Calc: GDAL_Calc.py -A Mask.tif -B CutBigImageToClip.tif --outfile=SmallerFile.tif --NoDataValue=0 --Calc="B*(A>0)" This query will populate 0 where Mask.tif <= 0 and BigImage where the Mask > 0. To do this both rasters must be the same cell size,...


8

Use gdaldem to add color to your tif file. gdaldem color-relief input.tif style.txt -alpha output.tif style.txt looks like this, each line having up to five values. First one is the raster's pixel value, next three are RGB values, and the fifth one is alpha, which is used for transparency - 255 means (by default) no transparency, 0 denotes fully transparent....


8

I would suggest trying to add the -et (error threshold) option with a lower thresholds than the default. If you try a lower -et threshold the horizontal artifact should disappear Or try changing the re-sampling gdalwarp -r mode As a reference gdalwarp


8

Using the GDAL Python bindings (GDAL >= 2.1), it should be: from osgeo import gdal input_raster = "path/to/rgb.tif" output_raster = "path/to/rgb_output_cut.tif" input_kml = "path/to/extent.kml" ds = gdal.Warp(output_raster, input_raster, format = 'GTiff', cutlineDSName = input_kml, cutlineLayer = '...


7

You can either leave the -s_srs empty or use +proj=sinu +R=6371007.181 +nadgrids=@null +wktext gdalwarp -of GTIFF -s_srs '+proj=sinu +R=6371007.181 +nadgrids=@null +wktext' -r cubic -t_srs '+proj=longlat +datum=WGS84 +no_defs' {inputfile} {outputfile} Note that inputfile must be the full subdataset (SDS) name, which you can get from gdalinfo. For example, ...


7

I think all you need to do is reproject it using: gdalwarp -co TILED=YES -co COMPRESS=DEFLATE -t_srs EPSG:3857 image.tif newImage.tif and then tile it: gdal2tiles.py newImage.tif If your file is very large it make take a while.


7

(Option 1) You could try the option in gdalwarp to ignore invalid polygons. This is invoked using "--config GDALWARP_IGNORE_BAD_CUTLINE YES": gdalwarp -of GTiff -co "COMPRESS=DEFLATE" -tr 30 30 -cutline VectorName -crop_to_cutline InputRaster OutputRaster --config GDALWARP_IGNORE_BAD_CUTLINE YES (Option 2) Since your polygon is a square, you could read the ...


7

It turns out you can actually use default parameter values calling Processing algorithms from PyQGIS using a different syntax. I didn't find it in the docs, but that's what we have GIS.SE for :D. Just call the algorithm providing parameters as a Python dictionary with keys being parameter names. A minimal example: import processing processing.runalg( "...


7

I had the same issue. Might be a bug in gdal configuration or something. To fix it, do this: sudo ln -s /usr/lib/grass72/lib/libgrass_* /usr/lib/


7

The best solution for you would be the VSIMEM filesystem which lets you save outputs of gdal utilities into a filesystem stored in memory. gdal.Warp('/vsimem/reprojected.tif', ds, srsSRS=in_proj, dstSRS=out_proj) You could also make a vrt: gdal.Warp('/vsimem/reprojected.vrt', ds, srsSRS=in_proj, dstSRS=out_proj) Once stored in the vsimem filesystem, the ...


6

gdal_merge.py loads all files into memory before processing them. therefore it is not able to process large files if your memory is small. see here


6

Here's a world file for you (as per Gene's suggestion). Cut and paste the following into a text editor and name the file the same as your full-sized map but give it the extension '.pgw' instead on '.png': 1.0 0.0 0.0 -1.0 0.0 3231.0 I've made the following assumptions: Your map has a resolution (pixel size) of 1m square (i.e. each pixel is a reasonable ...


6

For the Mercator projection, the extent can not reach North and South pole for mathematical reasons. The standard Google and Openstreetmap mercator projection is limited to 85.011° North and South to get a square map. See http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#X_and_Y for explanation. Using EPSG:3857, the extent of a map is -20037508,-...


Only top voted, non community-wiki answers of a minimum length are eligible