Spatial subset of geotiff through masking by shapefile?

I would like to create a mask (using poly2mask) to apply to a geotiff image, based on a single polygon in shapefile .shp format.

I have put together simple code below but evidently I have something wrong with syntax (maybe even just totally improper usage) of the worldtoIntrinsic part.

I have searched help avenues, but haven't cracked it.

Error message: (An issue with the structure of R, that I am returning by geotiffread?)

Undefined function 'worldToIntrinsic' for input arguments of type 'map.rasterref.GeographicPostingsReference'. Error in polygon_shapefileread_example (line 15) [ix, iy] = worldToIntrinsic(R,rx,ry);

Attached are the small geotiff and shapefile.

Many thanks for any help, Matt

 %Read in GeoTIFF
[I R] = geotiffread('geotiff_example.tif');
 % Read Region of Interest shapefile
roi = shaperead('shapefile_example.shp');
 % Remove trailing nan from shapefile
rx = roi.X(1:end-1);
ry = roi.Y(1:end-1);

 % convert to image coordinates
[ix, iy] = worldToIntrinsic(R,rx,ry);
 %Make the mask
mask = poly2mask(ix,iy,R.RasterSize(1),R.RasterSize(2));
 %Following line checks some 1's are generated in mask
maskcheck=sum(sum(mask));

Answer 1

[Chad Greene]

Does this work for you?

% Get image info: 
R = geotiffinfo('geotiff_example.tif'); 

% Get x,y locations of pixels: 
[x,y] = pixcenters(R); 

% Convert x,y arrays to grid: 
[X,Y] = meshgrid(x,y);

roi = shaperead('shapefile_example.shp');

% Remove trailing nan from shapefile
rx = roi.X(1:end-1);
ry = roi.Y(1:end-1);

mask = inpolygon(X,Y,rx,ry);

[mcb001]

Chad, this was an excellent piece of help, thanks a lot.

[Chad Greene]

Glad to hear it!

[Elias Berra]

The code works perfect for a single polygons, thanks for that. If my shapefile has multiple polygons, how could I create a mask containing the multiple polygons?

[Chad Greene]

Elias - You might try Kelly Kearney's inpolygons function.

[Sujai Banerji]

When I'm using the following code, I get this error:

% Remove trailing nan from shapefile rx = roi.X(1:end-1);

Error using end Incorrect cell or structure reference involving "end". Most likely cause is a reference to multiple elements of a cell or structure followed by additional subscript or structure references.

[Muhammad Habibie]

I try using them not using tif file but matrix file and follow your code, I get this error :

Error using parsePixMapInputs>parseOther (line 66) Function PIXCENTERS expected input number 1, R, to be either a 3-by-2 referencing matrix or a scalar map raster reference object. Instead its size was 9600x19200.

Error in parsePixMapInputs (line 22) [R, height, width, hasFlag] = parseOther(func_name, flagstr, varargin{:});

Error in pixcenters (line 51) = parsePixMapInputs('PIXCENTERS', 'makegrid', varargin{:});

Error in crop_matrix (line 74) [x,y] = pixcenters(I11);

Need suggestion.

[Reema Alhassan]

hello chad your solution works with me but I need to use it with large geotiff about 4 GB and when I run the code it crashes every time in meshgrid() or inpolygon() do you have any idea of how to deal with this ?

Thank you

[Veronica Bozzini]

Muhammad did you solve your problem? I get the same error.

[Abhishek Chakraborty]

Hi Chad. How to get the subset of the parent geotiff file ( R) ?

[André Luiz Lourenço]

@Abhishek Chakraborty to obtain a subset of Spatial Reference "R" you can use mapcrop funcion.

The function is explained here: https://www.mathworks.com/help/map/ref/mapcrop.html

[kinga kulesza]

for multiple polygons:

[X,Y]=meshgrid(lon,lat);
roi = shaperead('multiple_polygons.shp'); 
for i=1:numel(roi)
    rx = roi(i).X(1:end-1);
    ry = roi(i).Y(1:end-1);
    mask(:,:,i) = inpolygon(X,Y,rx,ry); 
end