Wie kann der Umfang der Rasterverarbeitung mit einer räumlichen Maske eingeschränkt werden?

Question

Ich versuche, die Rasterverarbeitung in MATLAB so zu beschränken, dass nur Bereiche innerhalb einer Shapefile-Grenze eingeschlossen werden, ähnlich wie ArcGIS Spatial Analyst-Funktionen eine verwenden. Hier finden Sie einige (reproduzierbar) Beispieldaten mit denen ich arbeite:

• A 4-band NAIP image (WARNING 169MB Download)
• A shapefile of study area boundaries (A Reißverschluss Shape-Datei auf File Dropper)

Hier ist ein MATLAB-Skript ich benutze NDVI zu berechnen:

file = 'C:\path\to\doi1m2011_41111h4nw_usda.tif'; 
[I R] = geotiffread(file); 
outputdir = 'C:\output\' 

% Calculate NDVI 
NIR = im2single(I(:,:,4)); 
red = im2single(I(:,:,1)); 

ndvi = (NIR - red) ./ (NIR + red); 
double(ndvi); 
imshow(ndvi,'DisplayRange',[-1 1]); 

% Stretch to 0 - 255 and convert to 8-bit unsigned integer 
ndvi = floor((ndvi + 1) * 128); % [-1 1] -> [0 256] 
ndvi(ndvi < 0) = 0;    % not really necessary, just in case & for symmetry 
ndvi(ndvi > 255) = 255;   % in case the original value was exactly 1 
ndvi = uint8(ndvi);    % change data type from double to uint8 

% Write NDVI to .tif file (optional) 
tiffdata = geotiffinfo(file); 
outfilename = [outputdir 'ndvi_' 'temp' '.tif']; 
geotiffwrite(outfilename, ndvi, R, 'GeoKeyDirectoryTag', tiffdata.GeoTIFFTags.GeoKeyDirectoryTag) 

Das folgende Bild zeigt, was ich möchte mit MATLAB erreichen. Für dieses Beispiel habe ich ArcGIS raster calculator (Float (Band4-Band1)/Float (Band4 + Band1)) verwendet, um den NDVI auf der rechten Seite zu erzeugen. Ich habe auch den Studienbereich Shapefile als mask in the environment settings angegeben.

Frage:

Wie kann ich die Rasterverarbeitung Ausmaß in MATLAB begrenzen eine Polygon Shape-Datei als räumliche Maske unter Verwendung der Ergebnisse in der Abbildung dargestellt zu replizieren?

Was ich habe erfolglos versucht:

roipoly und poly2mask, obwohl ich diese Funktionen nicht richtig anwenden kann, scheint (unter Berücksichtigung dieser räumlichen Daten sind), um die gewünschten Effekte zu erzeugen.

EDIT:

Ich habe versucht, die im Anschluss an die Shape-Datei auf eine Maske zu konvertieren, ohne Erfolg. Nicht sicher, wohin ich gehe falsch hier ...

s = 'C:\path\to\studyArea.shp' 

shp = shaperead(s) 
lat = [shp.X]; 
lon = [shp.Y]; 

x = shp.BoundingBox(2) - shp.BoundingBox(1) 
y = shp.BoundingBox(3) - shp.BoundingBox(1) 

x = poly2mask(lat,lon, x, y) 

---

Fehlermeldungen:

Error using poly2mask 
Expected input number 1, X, to be finite. 

Error in poly2mask (line 49) 
validateattributes(x,{'double'},{'real','vector','finite'},mfilename,'X',1); 

Error in createMask (line 13) 
x = poly2mask(lat,lon, x, y) 

Answer 1

Sie können die Region von Interesse gelesen von:

roi = shaperead('study_area_shapefile/studyArea.shp'); 

Hieb der hinteren NaN:

rx = roi.X(1:end-1); 
ry = roi.Y(1:end-1); 

Wenn Sie mehrere Polygone in Ihrem Shape-Datei haben, werden sie von NaNs getrennt und Sie haben um sie getrennt zu behandeln.

dann die worldToIntrinsic-Methode aus dem räumlichen Bezugs des sat-Bild verwenden, um die Polygon-Punkte in Bildkoordinaten zu konvertieren:

[ix, iy] = R.worldToIntrinsic(rx,ry); 

Dies setzt voraus, beide Koordinatensysteme sind die gleichen.

Dann können Sie gehen und Ihre Maske machen durch:

mask = poly2mask(ix,iy,R.RasterSize(1),R.RasterSize(2)); 

Sie die Maske auf Ihrem ursprünglichen mehrschichtige Bild verwenden können, bevor sie eine Berechnung machen von:

I(repmat(~mask,[1,1,4])) = nan; 

Oder es auf einem einzigen verwenden Schicht (zB rot) von:

red(~mask) = nan; 

Wenn die Regionen sehr klein sind, könnte es von Vorteil sein, (für Speicher und Rechenleistung), um ein maskiertes Bild in eine dünn besetzte Matrix umzuwandeln. Ich habe nicht versucht, wenn das einen Geschwindigkeitsunterschied macht.

red(~mask) = 0; 
sred = sparse(double(red)); 

Unfortunatly, spärlich Matrizes sind nur mit Doppel, so dass Ihr uint8 vor muss umgewandelt werden.

Im Allgemeinen sollten Sie die ROI aus dem Bild zuschneiden. Suchen Sie in den Objekten "roi" und "R" nach nützlichen Parametern und Methoden. Ich habe es hier nicht getan.

Schließlich meine Version des Skripts, mit einigen leichten anderen Veränderungen:

file = 'doi1m2011_41111h4nw_usda.tif'; 
[I R] = geotiffread(file); 
outputdir = ''; 

% Read Region of Interest 
roi = shaperead('study_area_shapefile/studyArea.shp'); 
% Remove trailing nan from shapefile 
rx = roi.X(1:end-1); 
ry = roi.Y(1:end-1); 
% convert to image coordinates 
[ix, iy] = R.worldToIntrinsic(rx,ry); 
% make the mask 
mask = poly2mask(ix,iy,R.RasterSize(1),R.RasterSize(2)); 
% mask sat-image 
I(repmat(~mask,[1,1,4])) = 0; 

% convert to sparse matrizes 
NIR = sparse(double(I(:,:,4))); 
red = sparse(double(I(:,:,1))); 
% Calculate NDVI 
ndvi = (NIR - red) ./ (NIR + red); 
% convert back to full matrizes 
ndvi = full(ndvi); 
imshow(ndvi,'DisplayRange',[-1 1]); 

% Stretch to 0 - 255 and convert to 8-bit unsigned integer 
ndvi = (ndvi + 1)/2 * 255; % [-1 1] -> [0 255] 
ndvi = uint8(ndvi);   % change and round data type from double to uint8 

% Write NDVI to .tif file (optional) 
tiffdata = geotiffinfo(file); 
outfilename = [outputdir 'ndvi_' 'temp' '.tif']; 
geotiffwrite(outfilename, ndvi, R, 'GeoKeyDirectoryTag', tiffdata.GeoTIFFTags.GeoKeyDirectoryTag); 
mapshow(outfilename); 

Answer 2

Es gibt drei Schritte hier, für die ich drei Funktionen erstellen:

1. Berechnen Sie das NDVI für das vollständige Eingabebild: ndvi = comp_ndvi(nir, red)
2. Berechnen Sie die m von der Shape-Datei fragen: mask = comp_mask(shape)
3. die NDVI Kombinieren und die Maske: output = combine_ndvi_mask(ndvi, mask)

Sie haben den Code für comp_ndvi() in Ihrer Frage.Der Code für combine_ndvi_mask() hängt davon ab, was Sie mit den maskierten Bereichen tun möchten. wenn man sie weiß machen wollen, könnte es wie folgt aussehen:

function output = combine_ndvi_mask(ndvi, mask) 
output = ndvi; 
output(~mask) = 255; 
end 

In comp_mask() wollen Sie poly2mask() verwenden, um die Polygon-Scheitelpunkte in der Rastermaske zu konvertieren. Um hier zu helfen, muss ich wissen, was Sie bereits haben. Haben Sie die Scheitelpunkte in MATLAB geladen? Was hast du mit Poly2mask versucht?