2

I need to clip my raster based on a polygonal shapefile, basically leave raster values inside the polygon shape (inside red perimeter in the figure below) as shown here:

enter image description here

Therefore, I use:

gdalwarp -cutline border200millas_polygon.shp -crop_to_cutline chl.2015.04.01_2015.04.08.aqua-color.tiff out.tiff

This should do the trick, but it doesn't. As result I get a black raster. After checking and comparing output values from both original and clipped raster I noticed that raster values of the clipped raster were much lower that those of the original raster (e.g., STATISTICS_STDDEV in out.tiff, this one being the clipped raster):

Files: chl.2015.04.01_2015.04.08.aqua-color.tiff                           |    Files: out.tiff
Size is 2606, 2553                                                         |    Size is 1632, 1913
Origin = (-90.463892341650677,0.065232508032915)                           |    Origin = (-84.673288999992508,-3.392115000001818)
Pixel Size = (0.008784683301344,-0.008784616065831)                        |    Pixel Size = (0.008786027376501,-0.008785333088489)
Upper Left  ( -90.4638923,   0.0652325) ( 90d27'50.01"W,  0d 3'54.84"N)    |    Upper Left  ( -84.6732890,  -3.3921150) ( 84d40'23.84"W,  3d23'31.61"S)
Lower Left  ( -90.4638923, -22.3618923) ( 90d27'50.01"W, 22d21'42.81"S)    |    Lower Left  ( -84.6732890, -20.1984572) ( 84d40'23.84"W, 20d11'54.45"S)
Upper Right ( -67.5710077,   0.0652325) ( 67d34'15.63"W,  0d 3'54.84"N)    |    Upper Right ( -70.3344923,  -3.3921150) ( 70d20' 4.17"W,  3d23'31.61"S)
Lower Right ( -67.5710077, -22.3618923) ( 67d34'15.63"W, 22d21'42.81"S)    |    Lower Right ( -70.3344923, -20.1984572) ( 70d20' 4.17"W, 20d11'54.45"S)
Center      ( -79.0174500, -11.1483299) ( 79d 1' 2.82"W, 11d 8'53.99"S)    |    Center      ( -77.5038907, -11.7952861) ( 77d30'14.01"W, 11d47'43.03"S)
Band 1 Block=2606x1 Type=Byte, ColorInterp=Red                             |    Band 1 Block=1632x1 Type=Byte, ColorInterp=Red
  Minimum=0.000, Maximum=255.000, Mean=160.832, StdDev=122.633             |      Minimum=0.000, Maximum=255.000, Mean=0.001, StdDev=0.479
    STATISTICS_MEAN=160.83171664774                                        |        STATISTICS_MEAN=0.00089845791949817
    STATISTICS_STDDEV=122.63291968596                                      |        STATISTICS_STDDEV=0.47865014597867
Band 2 Block=2606x1 Type=Byte, ColorInterp=Green                           |    Band 2 Block=1632x1 Type=Byte, ColorInterp=Green
  Minimum=0.000, Maximum=255.000, Mean=175.865, StdDev=108.243             |      Minimum=0.000, Maximum=255.000, Mean=0.001, StdDev=0.488
    STATISTICS_MEAN=175.86531037027                                        |        STATISTICS_MEAN=0.00096059725510696
    STATISTICS_STDDEV=108.24324530596                                      |        STATISTICS_STDDEV=0.48752743191956
Band 3 Block=2606x1 Type=Byte, ColorInterp=Blue                            |    Band 3 Block=1632x1 Type=Byte, ColorInterp=Blue
  Minimum=0.000, Maximum=255.000, Mean=252.883, StdDev=20.256              |      Minimum=0.000, Maximum=255.000, Mean=0.001, StdDev=0.540
    STATISTICS_MEAN=252.88313208935                                        |        STATISTICS_MEAN=0.0011434918975431
    STATISTICS_STDDEV=20.25592618995                                       |        STATISTICS_STDDEV=0.53998993166519

What would be the reason for this? Using -dstalpha just produces a blank output. I simply need to clip my original raster... nothing more than that. If I use gdal_rasterize instead, I can succesfully delete the polygonal shapefile from my original raster WITHOUT changing the raster original values, as explained here.

gdalinfo of the original raster below:

gdalinfo -mm -stats chl.2015.04.01_2015.04.08.aqua-color.tiff 
Driver: GTiff/GeoTIFF
Files: chl.2015.04.01_2015.04.08.aqua-color.tiff
Size is 2606, 2553
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-90.463892341650677,0.065232508032915)
Pixel Size = (0.008784683301344,-0.008784616065831)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -90.4638923,   0.0652325) ( 90d27'50.01"W,  0d 3'54.84"N)
Lower Left  ( -90.4638923, -22.3618923) ( 90d27'50.01"W, 22d21'42.81"S)
Upper Right ( -67.5710077,   0.0652325) ( 67d34'15.63"W,  0d 3'54.84"N)
Lower Right ( -67.5710077, -22.3618923) ( 67d34'15.63"W, 22d21'42.81"S)
Center      ( -79.0174500, -11.1483299) ( 79d 1' 2.82"W, 11d 8'53.99"S)
Band 1 Block=2606x1 Type=Byte, ColorInterp=Red
    Computed Min/Max=0.000,255.000
  Minimum=0.000, Maximum=255.000, Mean=160.832, StdDev=122.633
  Metadata:
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=160.83171664774
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=122.63291968596
Band 2 Block=2606x1 Type=Byte, ColorInterp=Green
    Computed Min/Max=0.000,255.000
  Minimum=0.000, Maximum=255.000, Mean=175.865, StdDev=108.243
  Metadata:
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=175.86531037027
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=108.24324530596
Band 3 Block=2606x1 Type=Byte, ColorInterp=Blue
    Computed Min/Max=0.000,255.000
  Minimum=0.000, Maximum=255.000, Mean=252.883, StdDev=20.256
  Metadata:
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=252.88313208935
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=20.25592618995
  • That could happen if your polygon is outside the image area but I guess that you have checked that. – user30184 Sep 17 '15 at 20:37
  • yeah thanks for suggesting that, I re-checked that, I will post an image of what I have as example. – Gery Sep 17 '15 at 20:49
  • Result is indeed odd. Could you prepare a zip that contains all the stuff needed for repeating your trial? – user30184 Sep 17 '15 at 21:40
  • perhaps it's the GDAL version, let me try 2.0 and let you know. – Gery Sep 18 '15 at 13:11
  • ok, I tried with GDAL 2.0.0 but I got the same result, if you get this working, could you share your succesful attempt? – Gery Sep 18 '15 at 18:39
3

I took one png file from the data_dir of GeoServer. Install GeoServer and you will have the same image to play with in directory:

\geoserver-2.7.1\data_dir\coverages\mosaic_sample\

Convert png into tiff and assign projection.

  gdal_translate -a_srs epsg:4326 global_mosaic_6.png  tiff.tiff
    Input file size is 50, 50
    0...10...20...30...40...50...60...70...80...90...100 - done.

Then use gdalwarp:

gdalwarp -cutline burn.shp -crop_to_cutline tiff.tiff cut.tif
Creating output file that is 16P x 15L.
Processing input file tiff.tiff.
Warning : the source raster dataset has a SRS, but the cutline features
not.  We assume that the cutline coordinates are expressed in the destination SRS. If not, cutline results may be incorrect.
0...10...20...30...40...50...60...70...80...90...100 - done.

This is the result:

enter image description here

Burn.shp contains one polygon:

POLYGON (( 10.375610236185853 39.384332051098106, 10.344426149120759 39.730475417520644, 10.87455562922735 39.92069834861771, 11.152094004106685 39.546489303836594, 10.92133175982499 39.16292503293594, 10.375610236185853 39.384332051098106 ))

It appeared later from the comments that the cutline shapefile "polygon border200millas_polygon.shp" contained three polygons and gdalwarp seems to accept only one geometry as a cutline. However, it can take a multipolygon which can be make from polygons for example with the QGIS function "dissolve".

  • AHA! perhaps it's because my polygon has 3 features? could you please add your ogrinfo -soto your answer? – Gery Sep 18 '15 at 18:54
  • also, you showed gdal_rasterize in your answer, how about gdalwarp? I also worked and produced the same output as with gdal_rasterize? – Gery Sep 18 '15 at 19:02
  • Ogrinfo is not interesting, it is just that one feature with no projection. But I guess you are right with the three polygons. Think, color white what is outside polygon A and then all that is outside polygon B -> everything colored. One 3-part MultiPolygon might work. – user30184 Sep 18 '15 at 19:06
  • I'm looking for a way to get a 1-feature polygon (instead of 3-features polygon), is that possible with GDAL? – Gery Sep 18 '15 at 19:08
  • it worked!! that was the issue, IT IS NEEDED A POLYGON WITH JUST 1 FEATURE!! please add the above statement to mark your answer as the answer. To get that, I used dissolve in QGIS, add also the way to dissolve multiparts polygon shape with GDAL. – Gery Sep 18 '15 at 19:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.