0

I am trying to trim a geotiff image by 500 pixels (500m) on all sides. The following is what I have tried so far, which reads a geotiff and describes various geotiff metrics such as the x-min and y-min of the bounding box. I have tried using imcrop(), although I get errors when trying to write the output to geotiff.

What is the best approach to trim the geotiff by 500 meters (or pixels) on all sides? 

%% Specify code and data locations
datadir = 'X:\temp\binary_data_temp\';     % set this to the directory where the data are stored; 1-band tiled images here
outputdir = 'X:\temp\temp_matlab_out2\';   % set this to the directory where files will be written

%%Load data
files = dir(fullfile(datadir, '*.tif'));
fileIndex = find(~[files.isdir]);

parfor ix = 1:length(fileIndex)         
    fileName = files(fileIndex(ix)).name;
    [Z, R]=geotiffread([datadir fileName]);

    %% Get geotiff metrics
    tiffdata = geotiffinfo([datadir fileName]);
    ncols = tiffdata.Width;
    nrows = tiffdata.Height;
    xllcorner = tiffdata.BoundingBox(1);
    yllcorner = tiffdata.BoundingBox(3);
    cellsize = tiffdata.UOMLengthInMeters;

    %% Crop image
    x500 = xllcorner + 500
    y500 = yllcorner + 500
    ncols500 = ncols - 1000
    nrows500 = nrows - 1000
    rectangle = [x500 y500 ncols nrows]

    clipped = imcrop(Z, rectangle)
    clipped2 = uint8(clipped)

    %% Write arrays to geotiffs
    [flout] = strread(fileName, '%s', 'delimiter','.')

    % Write clipped to .tif  
    outfilename = [outputdir 'clipped_geotiff_' char(flout(1)) '.tif'];  
    geotiffwrite(outfilename, clipped2, R, 'GeoKeyDirectoryTag', tiffdata.GeoTIFFTags.GeoKeyDirectoryTag) 

end

These are the error message that I get when running this:

>> clip_raster
Starting parallel pool (parpool) using the 'local' profile ... connected to 8 workers.
Error using geotiffwrite>validateImage (line 383)
Expected input number 2, A, to be nonempty.

Error in geotiffwrite>validateInputs (line 337)
[A, cmap] = validateImage(A, cmap);

Error in geotiffwrite (line 235)
[filename, A, cmap, R, Params] = validateInputs( ...

Error in clip_raster (line 16)
parfor ix = 1:length(fileIndex)         % Note this is the MATLAB parallel processing version of a for loop
Borealis
  • 8,044
  • 17
  • 64
  • 112
  • 1
    Why imcrop does not work. It does exactly that. – Ander Biguri Feb 05 '15 at 17:46
  • 1
    imcrop should work, remember that you need to adjust the header data to match when you write it out. – Raab70 Feb 05 '15 at 17:48
  • @Raab70 and Ander Biguri I have edited the post to include my implementation of the crop and the error messages. – Borealis Feb 05 '15 at 18:33
  • 1
    The error looks like it is during validation of the image header data. After you do your cropping and before you write out the image, update tiffdata.Width, height and bounding box to reflect the changes you made. Otherwise it's an invalid header. There may be other header attributes you need to update as well. – Raab70 Feb 05 '15 at 18:37
  • 1
    Take a look at the help documentation on [`geotiffwrite`](http://www.mathworks.com/help/map/ref/geotiffwrite.html), one of the examples is this exact problem (the last one). Notice that they create new metadata, they call it `R` and the new one is `subR` for the smaller image, before writing it out. – Raab70 Feb 05 '15 at 18:41

0 Answers0