The Stack Overflow podcast is back! Listen to an interview with our new CEO.
33

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.


29

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


21

You could explicitly set the output coordinate range using the target extent option to gdalwarp (ie. "-te -180 -90 180 90") but you can also use the CENTER_LONG configuration option to force rewrapping around a new central longitude. Something like this: gdalwarp -t_srs WGS84 ~/0_360.tif 180.tif -wo SOURCE_EXTRA=1000 \ --config CENTER_LONG 0 ...


20

While I don't know why GDAL provides this overlap in functionality, be sure to set the cache for gdalwarp to make it really fast: # assuming 3G of cache here: gdalwarp --config GDAL_CACHEMAX 3000 -wm 3000 $(list_of_tiffs) merged.tiff Be sure to not define more cache than having RAM on the machine.


20

Try to specify the nodata-value from your input raster and set it for the output as well. Furthermore add the option -crop_to_cutline to make exact crops. More about the options here. gdalwarp -srcnodata <in> -dstnodata <out> -crop_to_cutline -cutline INPUT.shp INPUT.tif OUTPUT.tif


16

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


16

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


12

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

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


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

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.


10

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


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


10

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


9

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'): ...


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


9

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


8

I would not recommend using the MODIS sinusoidal projection in analysis. It would be prudent to project your MODIS data to something a bit more tractable. You can request MODIS in a projected geographic (lat/long) coordinate system on the MODIS Golbal Subsets site . That said I have used this as my CRS for MODIS "+proj=sinu +R=6371007.181 +nadgrids=@null +...


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

Use galdem to add color to your tif file. gdaldem color-relief input.tif style.txt -alpha output.tif style.txt looks like this, with each line has five values. First one is value of raster, next three are RGB value, fifth one is alpha 3500 white 2500 235:220:175 50% 190 185 135 700 240 250 150 0 50 180 50 nv 0 0 0 0 By default, ...


8

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


7

I have a vague outline of how to do this, but there's plenty that you may have to understand. NetCDF files are complex general data containers so its not always clear how to get spatial data out of them. In this case, you can get the Soil_Moisture variable and that is just a 2d matrix with no coordinate reference. If you do image(A) you should see your soil ...


7

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


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


6

Basically you need to cut the raster into two parts and piece them back together with a new offset/scale. There's an example here of how to do that from [-180,180] to [0,360] with gdal_translate and the VRT driver: http://trac.osgeo.org/gdal/wiki/UserDocs/RasterProcTutorial Scan down to the "5 min tutorial" and the details are under "Virtual Files". It ...


6

Your source coordinate system is most likely not defined in the CVS file that GDAL searches for proj4 strings. It looks like you might be able to pass the source EPSG as 3031 (from spatialreference.org) Note that it looks like your input is in a local projection. Is this clipped from a larger raster? To explicitly define the source you could just provide ...


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

I can reproduce the error with another Natural Earth image http://naciscdn.org/naturalearth/10m/raster/NE2_LR_LC.zip and GDAL version 2.0-dev. Write a mail into gdal-dev mailing list. I believe that the issue is real and you can file next a ticket into GDAL bug tracker. I also verified that problem is not caused by the lanczos resampling but "-r average" ...


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