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:Wie kann der Umfang der Rasterverarbeitung mit einer räumlichen Maske eingeschränkt werden?
- 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)
Es hängt von der Verarbeitung. Welche Art von Funktion zum Beispiel? – chappjc
@chappjc Zwei Verarbeitungsschritte in meinem Workflow umfassen die Berechnung von NDVI und die Ausführung eines Bildfilters (imfilter()) zur Berechnung der Baumabdeckung in%. – Borealis