How do i extract 3 bands out of a 4 band image using NumpyArrayToRaster()?

Question

I need to extract 3 bands from a 4 band image. I am using a function called NumpyArrayToRaster() it accepts only upto 3 band images . how do i make it work for 4 band images ?

This is my code right now-

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')

The problem i am facing is with the conversion of numpy array back to raster image. The function NumpyArrayToRaster() accepts upto a 3 band image. I am currently working with a 4 band image and have to convert this to raster. So, i need suggestions on how to convert this to a 3 band image or any other function that does not have any band restrictions. — Avinash Paritala, Oct 25 '15 at 18:10
http://pastie.org/10507013# is the link to the code i have tried — Avinash Paritala, Oct 25 '15 at 18:11
are the target 3 bands with different wavelength(range) then in the 4 band image ? If yes then you need to reintegrate the bands back .... see [Multi-Band Image raster to RGB](http://stackoverflow.com/a/29575362/2521214) — Spektre, Oct 26 '15 at 08:44

score 0 · Answer 1 · answered Feb 08 '18 at 12:26

As reading your code, I realized that you need to calculate NDVI values. Maybe you con slightly modify your approach and instead of using function NumpyArrayToRaster() you can use a simpler approach?

Here I offer my code, which: - reads the individual stacks of Sentinel data from multiple dates (multiband composites) - identify bands needed for NDVI calculation using rasterName/Band_X - calculate NDVI from bands - save NDVIs for each date to output folder, keep its name and date

In Landsat, those bands for NDVI are Band_4 and Band_3

My code:

# 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"

How do i extract 3 bands out of a 4 band image using NumpyArrayToRaster()?

1 Answers1