Trimimg a geotiff by a shapefile

Hi,

I am new to mapping toolbox and here the thing I need help: I have a large geotiff file and I like to extract a small area of it (still in a geotiff format) based on a shape file I have access to (I tried 'maptrims' but it requires the 'RasterInterpretation' field of the 'spatial referencing object' to be 'cells' which in my case it is 'postings'). Any idea on how can I achieve this?

Thanks,

Answer 1

[Chad Greene]

You can use geoimread to read only the data surrounding your shapefile. Get the xdata and ydata from the shapefile, then use those coordinates with geoimread.

[Mohammad Gohardoust]

Thanks Chad for your response. Finally I did it by 'maptrims' and then employing 'vec2mtx' but using 'geoimread' sounds more wise to not read the whole image data.

[Reema Alhassan]

Hello Mohammad could you share your code please if you still have it

thank you .

[Mohammad Gohardoust]

Hello Reema,

Assuming the map and shape files have geographic coordinates, one approach would be:

shape_name = 'shape_name.shp';
map_name = 'map_name.tif';
roi = shaperead(shape_name);
[A, R] = geotiffread(map_name);

shape_lon = roi.X;
shape_lat = roi.Y;

[A1,~] = vec2mtx(shape_lat, shape_lon, A, R, 'filled');
[row, col] = find(A1 ~= 2);
[lat, lon] = pix2latlon(R, row, col);
latlim = [min(lat)-R.CellExtentInLatitude/2, max(lat)+R.CellExtentInLatitude/2];
lonlim = [min(lon)-R.CellExtentInLongitude/2, max(lon)+R.CellExtentInLongitude/2];

A2 = A(min(row):max(row), min(col):max(col));
R2 = georefcells(latlim,lonlim,size(A2));
geotiffwrite('out_mapname.tif', A2, R2);

[Reema Alhassan]

Thanks Mohammad !! but do you know if I have a shape file with many polygons how can I clip the image for each polygon ? I mean the x and y coordinates from the shape file for the hole shape but now I need to specify the coordinates if you have an idea of doing this I will be very thankful

Many thanks,

[Mohammad Gohardoust]

I'm not sure whether or not you found the answer yet, but if I got the question correctly, multiple polygons in shape files are (usually I believe) separated by a 'NaN'. Using this fact you can put a for loop after reading the shape file data, with nnz(isnan(roi.X))+1 repetitions.