How do I save 6bands layer stacked data in ENVI format?

I have layer stacked 6 bands of sentinel2 data which is in jp2 format after resizing green,red and nir to 20m resolution. I want to save this 6bands layer stacked data into ENVI format. I request you all to please suggest me how to do it.

I would appreciate your kind help.

Devendra

Answer 1

[Vinayak Choyyan]

Hi Devendra,

Please check out the below sample code to read, layer stack, and then write Sentinel 2 images as ENVI files.

%Read 6 bands
B02=imread("R10m\tile_id_B02_10m.jp2");
B03=imread("R10m\tile_id_B03_10m.jp2");
B04=imread("R10m\tile_id_B04_10m.jp2");
B08=imread("R10m\tile_id_B08_10m.jp2");
B11=imread("R20m\tile_id_B11_20m.jp2");
B12=imread("R20m\tile_id_B12_20m.jp2");
%resize desired bands
bandSize20m=5490;
B02=imresize(B02,[bandSize20m,bandSize20m]);
B03=imresize(B03,[bandSize20m,bandSize20m]);
B04=imresize(B04,[bandSize20m,bandSize20m]);
B08=imresize(B08,[bandSize20m,bandSize20m]);
%layer stack the bands
cube=zeros(bandSize20m,bandSize20m,6,"uint16");%6 bands
cube(:,:,1)=B02;
cube(:,:,2)=B03;
cube(:,:,3)=B04;
cube(:,:,4)=B08;
cube(:,:,5)=B11;
cube(:,:,6)=B12;
%set wavelength of each band
wavelength=[490 560 665 842 1610 2190];
%create metadata struct
metadata.WavelengthUnits="nm";
%create hypercube object
hcube=hypercube(cube,wavelength,metadata);
%write as ENVI file
enviwrite(hcube,"location\filename","DataType","uint16")

In the above demo code, I have only set one metadata as an example. You can set more metadata fields by adding to the stuct named 'metadata' in above code and all of those will apprear in the .hdr file generated by enviwrite.

You can read more about parameterised hypercube object creation here and also check out enviwrite function here.

Hope this helps. Regards.

[Vinayak Choyyan]

Additionally, I advice adding keywords like 'hyperspectral' or 'multispectral' to tags. This would ensure prompt notification and response to your inquiry.

[gauri]

Thank you very much for your kind help. I want to add map and projection information in ENVI header file as follows;

map info = {UTM, 1.000, 1.000, 712580.000, 3216500.000, 2.0000000000e+001, 2.0000000000e+001, 43, North, WGS-84, units=Meters}

coordinate system string = {PROJCS["WGS_1984_UTM_Zone_43N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",75.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]}

Please suggest me how to add this information in ENVI header file.

Devendra

[Vinayak Choyyan]

For MapInfo add the following to metadata struct

metadata.MapInfo.ProjType="UTM";
metadata.MapInfo.PixelTiePoints=[1 1];
metadata.MapInfo.MapTiePoints=[712580.000, 3216500.000];
metadata.MapInfo.PixelSize=[20 20];
metadata.MapInfo.ProjZone=43;
metadata.MapInfo.NorthOrSouth="North";
metadata.MapInfo.Datum="WGS-84";
metadata.MapInfo.Units="Meters";

and for coordinate system string add the following to metadata struct

%add the entire coordinate system string data as a character array
metadata.CoordinateSystemString='{PROJCS["WGS_1984_UTM_Zone_43N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",75.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]}';

This should give you the above to info in ENVI .hdr file when using enviwrite().

[gauri]

Thanks a lot for your kind help. However, it has not written ENVI header file as it should appear in ENVI header file. I am pasting the generated header file as follows;

ENVI

bands = 7

data type = 12

interleave = bsq

header offset = 0

coordinate system string = {{PROJCS["WGS_1984_UTM_Zone_43N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",75.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]}}

wavelength units = Nanometers

byte order = 0

wavelength = {560, 665, 705, 740, 783, 842, 865}

lines = 5490

samples = 5490

map info = {UTM, 1, 1, 712580, 3216500, 20, 20, 43, North, WGS-84, Meters}

Further, I have stored the names of files which are layer stacked in cube as band_names. How do I write band_names in ENVI header file also?.

[Vinayak Choyyan]

This is a valid ENVI header file and you should be able to import this to ENVI. I am not sure what you meant by "as it should appear in ENVI header file".

For band names,

metadata.BandNames=["B02","B03","B04"];

Any other info that you would like to add to header file can be added using the metadata struct as I have mentioned, that is, by adding it to the struct.

[gauri]

Thank you once again for your support. It worked well without BandNames. Could you please suggest me how to add BandNames in metadata? I tried as per your suggestion but it did not work.

I appreciate your kind cooperation and help.

[Vinayak Choyyan]

Please add more details as to what did not work. It will help me understand at which step you are facing an issue. Please try this for adding band names

metadata.BandNames=["B02","B03","B04","B08","B11","B12"].';

[gauri]

If I add metadata.BandNames=["B02","B03","B04","B08","B11","B12"].'; than following errors occured;

Error using uint8

Conversion to uint8 from cell is not possible.

Error in hypercube>array2Stream (line 2078)

stream = uint8(array);

Error in hypercube>getHeaderStream (line 2047)

data = array2Stream(thisField);

Error in hypercube>enviheaderwrite (line 1979)

stream = getHeaderStream(header);

Error in hypercube/enviwrite (line 971)

enviheaderwrite(obj.Metadata, hdrFile);

Error in layer_stack (line 56)

enviwrite(hcube,fname,"DataType","uint16")

Please suggest me how to fix it.

[gauri]

If I add metadata.BandNames=["B02","B03","B04","B08","B11","B12"].'; than following errors occured;

Error using uint8

Conversion to uint8 from cell is not possible.

Error in hypercube>array2Stream (line 2078)

stream = uint8(array);

Error in hypercube>getHeaderStream (line 2047)

data = array2Stream(thisField);

Error in hypercube>enviheaderwrite (line 1979)

stream = getHeaderStream(header);

Error in hypercube/enviwrite (line 971)

enviheaderwrite(obj.Metadata, hdrFile);

Error in layer_stack (line 56)

enviwrite(hcube,fname,"DataType","uint16")

Please suggest me how to fix it.

[Vinayak Choyyan]

This is not expected, and I do not get this error. I had tried the code before I shared it with you. Please share your code where this error is reproduceable, and I will have a look.

[gauri]

I am attaching the code and request you to please have a look on it and suggest me how to add band names also in ENVI header file.

I also want to extract the data using shape file from this seven band stacked layer file which is a seperate question on following link;

https://in.mathworks.com/matlabcentral/answers/2116911-how-do-i-extract-the-multiband-data-using-shape-file?s_tid=mlc_ans_email_ques

I request you to please have a look on it and help me.

[Vinayak Choyyan]

At line number 41, I see that you have passed the data as a row vector like this:

metadata.BandNames=  ["B03","B04","B05","B06","B07","B08","B8A"];

Please try passing it as a column vector like this as mentioned in my previous comments:

metadata.BandNames=  ["B03","B04","B05","B06","B07","B08","B8A"].';

Here I am using transpose operator. I have linked the documentation.

Additionally, please go through this example to see how to use a shapefile to crop and mask certain regions of your image data.