How to match GeoTiff file with Shapefile?

I have a Shapefile with a bunch of fields (Geometry, BoundingBox, Lon, Lat, etc.) where the geometry is in Polygon form. So I have a list of polygonal regions which I want to extract from my tiff files but I am having trouble doing this. I've read all the documentation that I could find on shapefiles and geotiff files but nothing has been able to help me answer this question.

I currently have the following code

[A, R] = geotiffread(...);

S = shaperead(..., 'UseGeoCoords', true, 'Selector', { ... });

So I have everything loaded into MatLab structs/matrices but I just can't figure out how to match up the two structures. How can I extract the data from A where the location is contained within the region(s) specified by S ?

Answer 1

[Chad Greene]

Two possible solutions:

Given some polygon described by the arrays lat and lon, you could use Aslak Grinsted's geoimread to load a rectangular region of the image that bounds the lat/lon array:

[A,x,y] = geoimread('filename.tif',lat,lon);

The above returns x/y coordinates of each pixel. Convert lat & lon to x/y coordinates with projfwd, then use inpolygon to get all the geotiff data within the polygon.

A different solution: If you just need to interpolate along a lat/lon array, or in some lat/lon grid, use geotiffinterp. Syntax would be

Ai = geotiffinterp('filename.tif',lat,lon);

[Ryan McKenna]

I tried the first function and got an error message:

No image coordinates can be found inside the specified limit

I also tried the second function but that doesn't work either. Any ideas what could be causing this?

[Chad Greene]

I'm guessing that your tiff file either does not have georeferencing metadata written into it, or it's there but geoimread and geotiffinterp can't find it. Or, perhaps your geotiff is properly referenced to some part of Europe while your lat/lon arrays outline Fiji. Yet another possibility is that your lon array contains longitudes of 0 to 360 while your tiff is referenced to -180 to 180 or vice versa.

[Elias Berra]

The tool works great. However, if you want extract information only from within a polygon (e.g. circle), that is not possible, since it consider the bounding box of the polygon instead the actual polygon limits.

[Chad Greene]

Elias, you'll have to use inpolygon for that.