Как извлечь 3 полосы из 4-полосного изображения с помощью NumpyArrayToRaster()?

Question

Мне нужно извлечь 3 полосы из 4-х полосного изображения. Я использую функцию NumpyArrayToRaster(), которая принимает только до 3-х полосных изображений. как мне заставить это работать для 4-полосных изображений?

Это мой код прямо сейчас

import arcpy
arcpy.CheckOutExtension("Spatial")
from PIL import Image

import matplotlib.pyplot as plt
import numpy as np
from skimage import io
from skimage.segmentation import quickshift

arcpy.env.overwriteOutput = True

# The input 4-band NAIP image
img = r'C:\Users\Alekhya\Desktop\Krishna\NRSC\processed image\newclip4\clip4.tif'


# Convert image to numpy array
imgarr = io.imread(img)
print imgarr
print imgarr.shape
print imgarr.dtype

# Run the quick shift segmentation
segments = quickshift(imgarr, kernel_size=3, convert2lab=False, max_dist=6, ratio=0.5)
print("Quickshift number of segments: %d" % len(np.unique(segments)))


# View the segments via Python
plt.imshow(segments)
print segments
print segments.shape
print type(segments)
print segments.dtype


# Get raster metrics for coordinate info
imgRaster = arcpy.sa.Raster(img)

# Lower left coordinate of block (in map units)
mx = imgRaster.extent.XMin
my = imgRaster.extent.YMin
sr = imgRaster.spatialReference

'''
# Note the use of arcpy to convert numpy array to raster
seg = arcpy.NumPyArrayToRaster(segments, arcpy.Point(mx,my), imgRaster.meanCellWidth, imgRaster.meanCellHeight)

outRaster = r'C:\Users\Alekhya\Desktop\Krishna\NRSC\processed image\newclip4\segments_clip4.tif'
seg_temp = seg.save(outRaster)
arcpy.DefineProjection_management(outRaster, sr)
'''

# Calculate NDVI from bands 4 and 3
b4 = arcpy.sa.Raster(r'C:\Users\Alekhya\Desktop\Krishna\NRSC\processed image\newclip4\clip4.tif\Band_4')
b3 = arcpy.sa.Raster(r'C:\Users\Alekhya\Desktop\Krishna\NRSC\processed image\newclip4\clip4.tif\Band_3')
ndvi = arcpy.sa.Float(b4-b3) / arcpy.sa.Float(b4+b3)
print ndvi

# Extract NDVI values based on image object boundaries
zones = arcpy.sa.ZonalStatistics(segments, "VALUE", ndvi, "MEAN")
zones.save(r'C:\Users\Alekhya\Desktop\Krishna\NRSC\processed image\newclip4\zones_clip4.tif')

# Classify the segments based on NDVI values
binary = arcpy.sa.Con(zones < 20, 1, 0)
binary.save(r'C:\Users\Alekhya\Desktop\Krishna\NRSC\processed image\newclip4\classified_clip4.tif')

Answer 1

Читая ваш код, я понял, что вам нужно рассчитать значения NDVI. Может быть, вы немного измените свой подход и вместо использования функции NumpyArrayToRaster(Вы можете использовать более простой подход?

Здесь я предлагаю свой код, который: - считывает отдельные стеки данных Sentinel из нескольких дат (многополосные композиты) - идентифицирует полосы, необходимые для расчета NDVI, используя rasterName/Band_X- рассчитать NDVI из диапазонов - сохранить NDVI для каждой даты в выходной папке, сохранить ее имя и дату

В Landsat эти полосы для NDVI - Band_4 и Band_3

Мой код:

# calculate NDVI from sentinel data
# list raster and calculate NDVI per each raster individually

# import modules
import arcpy, string

# import environmental settings
from arcpy import env
from arcpy.sa import *

# check out spatial extension
arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True

# add workspace
env.workspace = "C:/Users/input"

# List rasters
rasters = arcpy.ListRasters("*", "TIF")

# Output directory
outWd = "C:/Users/output/ndvi"

# calculate ndvi for every sentinel raster
for raster in rasters:
    # define ndvi outputs
    outNDVI = outWd + "/"+ raster.replace(".tif", "_ndvi.tif")
    print "outNDVI is " + outNDVI

    # specify inputs for ndvi and final output
    # NDVI takes NIR and Red, which are in Sentinel  Band 4 and Band 8 
    Red = raster + '\Band_4'
    NIR = raster + '\Band_8'

    # Create Numerator and Denominator rasters as variables, and
    # NDVI as output
    Num = arcpy.sa.Float(Raster(NIR) - Raster(Red))
    Denom = arcpy.sa.Float(Raster(NIR) + Raster(Red))
    NDVI = arcpy.sa.Divide(Num, Denom)
    print "NDVI calculating"

    # save results output
    NDVI.save(outNDVI)
    print "NDVI saved"