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.
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.
$ gdalwarp -tap -tr 30 30 -t_srs "etc..." -of ...
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.
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
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 ...
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 ...
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.
Create a reference to compare to:
gdalwarp -t_srs EPSG:4326 byte.tif warp_ref.tif
Run the ...
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.
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).
>cs2cs +init=EPSG:4326 +to +init=EPSG:3857
556597.45 5311971.85 0.00
1669792.36 6106854.83 0.00
the command should be something like:
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
to create each output raster.
Example of script:
# "shp" - folder for shapefiles
# "outputraster" - folder for output rasters
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 (...
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 ...
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):
for raster in findRasters(INPUT_FOLDER, '*.tif'):
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
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.
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....
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
50% 190 185 135
700 240 250 150
0 50 180 50
nv 0 0 0 0
By default, ...
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
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 ...
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,...
I think all you need to do is reproject it using:
gdalwarp -co TILED=YES -co COMPRESS=DEFLATE -t_srs EPSG:3857 newImage.tif image.tif
and then tile it:
If your file is very large it make take a while.
(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 ...
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:
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 ...