mapgen script with MODIS L2
Posted: Mon Jan 30, 2023 3:54 pm America/New_York
Hello,
I have been using the gpt command to process MODIS and VIIRS L2 images of chl and sst. I am using SeaDAS 8.3, with OCSSW installed, and Ubuntu 22.1. In the .xml (attached at the end of the post) in summary what I do is:
read(chl image) --> BandMaths(to create chl and flags bands) -->
Collocate(to bring all the bands together) --> Reproject --> Subset --> Write(as netcdf-CF)
read(sst image) --> BandMaths(to create chl and flags bands) -->
In another post, gnwiii recommended me that the correct treatment for L2 images is to use the mapgen function so as not to generate errors on the edges of the image. I have been reading about the mapgen function (https://seadas.gsfc.nasa.gov/help-8.3.0/processors/ProcessMapgen.html), but I was not able to execute it. When I executed the command: (I'm using Python 3.7)
python path/SeaDAS/ocssw/bin/mapgen --ifile=path/A2020003181000.L2_LAC_OC.x.nc --ofile=path/mapgenTest_A2020003181000.L2_LAC_OC.x.nc --product=chlor_a,l2_flags --resolution=1km --oformat=netcdf4 --projection=mercator --verbose
I get the error:
###### l2bin ######
Running: /home/fernando/SeaDAS/ocssw/bin/l2bin flaguse=ATMFAIL,CLDICE,BOWTIEDEL l3bprod=chlor_a,l2_flags ifile=/home/fernando/Escritorio/Pruebas_gpt/chl/A2020003181000.L2_LAC_OC.x.nc ofile=mapgen_tmp.l2bin.nc area_weighting=1 prodtype=regional resolution=H
Exiting: /home/fernando/SeaDAS/ocssw/bin/l2bin flaguse=ATMFAIL,CLDICE,BOWTIEDEL l3bprod=chlor_a,l2_flags ifile=/home/fernando/Escritorio/Pruebas_gpt/chl/A2020003181000.L2_LAC_OC.x.nc ofile=mapgen_tmp.l2bin.nc area_weighting=1 prodtype=regional resolution=H returned a status 1
I need output files in netcdf format (better if netcdf-CF, that's why I'm using gpt with Snap), what kind of netcdf is the output file, CF or BEAM-DIMAP?
My goal is to get the L2 image with a regular grid (no matter the projection), the bandmath steps like the subset can be skipped. I'm looking for the regular grid so I can work the images comfortably in python. Is the procedure I'm doing to process the images using the .xml correct?
Do you know any other way to get a regular grid from a L2 image?
thank you very much for your help
------------------------------------------------------
*file.xml*
<graph id="firstGraph">
<version>1.0</version>
<node id="Read_chl">
<operator>Read</operator>
<sources/>
<parameters class="com.bc.ceres.binding.dom.XppDomElement">
<file>${filechl}</file>
</parameters>
</node>
<node id="BandMaths_chl">
<operator>BandMaths</operator>
<sources>
<sourceProduct refid="Read_chl"/>
</sources>
<parameters>
<targetBands>
<targetBand>
<name>cldice_chl</name>
<expression>l2_flags.CLDICE</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>chlfail_chl</name>
<expression>l2_flags.CHLFAIL</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>land_chl</name>
<expression>l2_flags.LAND</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>bowtiedel_chl</name>
<expression>l2_flags.BOWTIEDEL</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>straylight_chl</name>
<expression>l2_flags.STRAYLIGHT</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>orig_chl</name>
<expression>chlor_a</expression>
<description></description>
<type>float32</type>
<validExpression></validExpression>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>chl</name>
<expression>if chlor_a > 8 or chlor_a < 0.06 then NaN else chlor_a</expression>
<description></description>
<type>float32</type>
<validExpression></validExpression>
<noDataValue>NaN</noDataValue>
</targetBand>
</targetBands>
</parameters>
</node>
<node id="Read_sst">
<operator>Read</operator>
<sources/>
<parameters class="com.bc.ceres.binding.dom.XppDomElement">
<file>${filesst}</file>
</parameters>
</node>
<node id="BandMaths_sst">
<operator>BandMaths</operator>
<sources>
<sourceProduct refid="Read_sst"/>
</sources>
<parameters>
<targetBands>
<targetBand>
<name>cldice_sst</name>
<expression>l2_flags.CLDICE</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>bad_sst</name>
<expression>qual_sst</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>land_sst</name>
<expression>l2_flags.LAND</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>bowtiedel_sst</name>
<expression>l2_flags.BOWTIEDEL</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>straylight_sst</name>
<expression>l2_flags.STRAYLIGHT</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>orig_sst</name>
<expression>sst</expression>
<description></description>
<type>float32</type>
<validExpression></validExpression>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>sst</name>
<expression>if sst > 30 or sst < 0 then NaN else sst</expression>
<description></description>
<type>float32</type>
<validExpression></validExpression>
<noDataValue>NaN</noDataValue>
</targetBand>
</targetBands>
</parameters>
</node>
<node id="CollocateNode">
<operator>Collocate</operator>
<sources>
<sourceProduct refid="BandMaths_chl"/>>
<sourceProduct.1 refid="BandMaths_sst"/>>
</sources>
<parameters class="com.bc.ceres.binding.dom.XppDomElement">
<targetProductType>COLLOCATED</targetProductType>
<resamplingType>NEAREST_NEIGHBOUR</resamplingType>
</parameters>
</node>
<node id="Reproject">
<operator>Reproject</operator>
<sources>
<sourceProduct refid="CollocateNode"/>
</sources>
<parameters class="com.bc.ceres.binding.dom.XppDomElement">
<wktFile/>
<crs>PROJCS["Mercator_1SP / World Geodetic System 1984",
GEOGCS["World Geodetic System 1984",
DATUM["World Geodetic System 1984",
SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
UNIT["degree", 0.017453292519943295],
AXIS["Geodetic longitude", EAST],
AXIS["Geodetic latitude", NORTH]],
PROJECTION["Mercator_1SP"],
PARAMETER["latitude_of_origin", 0.0],
PARAMETER["central_meridian", 0.0],
PARAMETER["scale_factor", 1.0],
PARAMETER["false_easting", 0.0],
PARAMETER["false_northing", 0.0],
UNIT["m", 1.0],
AXIS["Easting", EAST],
AXIS["Northing", NORTH]]</crs>
<resampling>Nearest</resampling>
<referencePixelX/>
<referencePixelY/>
<easting/>
<northing/>
<orientation/>
<pixelSizeX/>
<pixelSizeY/>
<width/>
<height/>
<tileSizeX/>
<tileSizeY/>
<orthorectify>false</orthorectify>
<elevationModelName/>
<noDataValue>NaN</noDataValue>
<includeTiePointGrids>true</includeTiePointGrids>
<addDeltaBands>false</addDeltaBands>
</parameters>
</node>
<node id="Subset">
<operator>Subset</operator>
<sources>
<sourceProduct refid="Reproject"/>
</sources>
<parameters class="com.bc.ceres.binding.dom.XppDomElement">
<sourceBands/>
<region>0,0,0,0</region>
<referenceBand/>
<geoRegion>POLYGON ((${east} ${north}, ${west} ${north}, ${west} ${south}, ${east} ${south}, ${east} ${north}))</geoRegion>
<subSamplingX>1</subSamplingX>
<subSamplingY>1</subSamplingY>
<fullSwath>false</fullSwath>
<tiePointGridNames/>
<copyMetadata>true</copyMetadata>
</parameters>
</node>
<node id="Write">
<operator>Write</operator>
<sources>
<source>Subset</source>
</sources>
<parameters>
<file>${ofile}</file>
<formatName>NetCDF4-CF</formatName>
</parameters>
</node>
</graph>
------------------------------------------------------
*file.par*
east=-51
west=-63
north=-35
south=-42
referenceComponentPattern=${ORIGINAL_NAME}_chl
dependentComponentPattern=${ORIGINAL_NAME}_sst
I have been using the gpt command to process MODIS and VIIRS L2 images of chl and sst. I am using SeaDAS 8.3, with OCSSW installed, and Ubuntu 22.1. In the .xml (attached at the end of the post) in summary what I do is:
read(chl image) --> BandMaths(to create chl and flags bands) -->
Collocate(to bring all the bands together) --> Reproject --> Subset --> Write(as netcdf-CF)
read(sst image) --> BandMaths(to create chl and flags bands) -->
In another post, gnwiii recommended me that the correct treatment for L2 images is to use the mapgen function so as not to generate errors on the edges of the image. I have been reading about the mapgen function (https://seadas.gsfc.nasa.gov/help-8.3.0/processors/ProcessMapgen.html), but I was not able to execute it. When I executed the command: (I'm using Python 3.7)
python path/SeaDAS/ocssw/bin/mapgen --ifile=path/A2020003181000.L2_LAC_OC.x.nc --ofile=path/mapgenTest_A2020003181000.L2_LAC_OC.x.nc --product=chlor_a,l2_flags --resolution=1km --oformat=netcdf4 --projection=mercator --verbose
I get the error:
###### l2bin ######
Running: /home/fernando/SeaDAS/ocssw/bin/l2bin flaguse=ATMFAIL,CLDICE,BOWTIEDEL l3bprod=chlor_a,l2_flags ifile=/home/fernando/Escritorio/Pruebas_gpt/chl/A2020003181000.L2_LAC_OC.x.nc ofile=mapgen_tmp.l2bin.nc area_weighting=1 prodtype=regional resolution=H
Exiting: /home/fernando/SeaDAS/ocssw/bin/l2bin flaguse=ATMFAIL,CLDICE,BOWTIEDEL l3bprod=chlor_a,l2_flags ifile=/home/fernando/Escritorio/Pruebas_gpt/chl/A2020003181000.L2_LAC_OC.x.nc ofile=mapgen_tmp.l2bin.nc area_weighting=1 prodtype=regional resolution=H returned a status 1
I need output files in netcdf format (better if netcdf-CF, that's why I'm using gpt with Snap), what kind of netcdf is the output file, CF or BEAM-DIMAP?
My goal is to get the L2 image with a regular grid (no matter the projection), the bandmath steps like the subset can be skipped. I'm looking for the regular grid so I can work the images comfortably in python. Is the procedure I'm doing to process the images using the .xml correct?
Do you know any other way to get a regular grid from a L2 image?
thank you very much for your help
------------------------------------------------------
*file.xml*
<graph id="firstGraph">
<version>1.0</version>
<node id="Read_chl">
<operator>Read</operator>
<sources/>
<parameters class="com.bc.ceres.binding.dom.XppDomElement">
<file>${filechl}</file>
</parameters>
</node>
<node id="BandMaths_chl">
<operator>BandMaths</operator>
<sources>
<sourceProduct refid="Read_chl"/>
</sources>
<parameters>
<targetBands>
<targetBand>
<name>cldice_chl</name>
<expression>l2_flags.CLDICE</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>chlfail_chl</name>
<expression>l2_flags.CHLFAIL</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>land_chl</name>
<expression>l2_flags.LAND</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>bowtiedel_chl</name>
<expression>l2_flags.BOWTIEDEL</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>straylight_chl</name>
<expression>l2_flags.STRAYLIGHT</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>orig_chl</name>
<expression>chlor_a</expression>
<description></description>
<type>float32</type>
<validExpression></validExpression>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>chl</name>
<expression>if chlor_a > 8 or chlor_a < 0.06 then NaN else chlor_a</expression>
<description></description>
<type>float32</type>
<validExpression></validExpression>
<noDataValue>NaN</noDataValue>
</targetBand>
</targetBands>
</parameters>
</node>
<node id="Read_sst">
<operator>Read</operator>
<sources/>
<parameters class="com.bc.ceres.binding.dom.XppDomElement">
<file>${filesst}</file>
</parameters>
</node>
<node id="BandMaths_sst">
<operator>BandMaths</operator>
<sources>
<sourceProduct refid="Read_sst"/>
</sources>
<parameters>
<targetBands>
<targetBand>
<name>cldice_sst</name>
<expression>l2_flags.CLDICE</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>bad_sst</name>
<expression>qual_sst</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>land_sst</name>
<expression>l2_flags.LAND</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>bowtiedel_sst</name>
<expression>l2_flags.BOWTIEDEL</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>straylight_sst</name>
<expression>l2_flags.STRAYLIGHT</expression>
<description></description>
<type>int32</type>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>orig_sst</name>
<expression>sst</expression>
<description></description>
<type>float32</type>
<validExpression></validExpression>
<noDataValue>NaN</noDataValue>
</targetBand>
<targetBand>
<name>sst</name>
<expression>if sst > 30 or sst < 0 then NaN else sst</expression>
<description></description>
<type>float32</type>
<validExpression></validExpression>
<noDataValue>NaN</noDataValue>
</targetBand>
</targetBands>
</parameters>
</node>
<node id="CollocateNode">
<operator>Collocate</operator>
<sources>
<sourceProduct refid="BandMaths_chl"/>>
<sourceProduct.1 refid="BandMaths_sst"/>>
</sources>
<parameters class="com.bc.ceres.binding.dom.XppDomElement">
<targetProductType>COLLOCATED</targetProductType>
<resamplingType>NEAREST_NEIGHBOUR</resamplingType>
</parameters>
</node>
<node id="Reproject">
<operator>Reproject</operator>
<sources>
<sourceProduct refid="CollocateNode"/>
</sources>
<parameters class="com.bc.ceres.binding.dom.XppDomElement">
<wktFile/>
<crs>PROJCS["Mercator_1SP / World Geodetic System 1984",
GEOGCS["World Geodetic System 1984",
DATUM["World Geodetic System 1984",
SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
UNIT["degree", 0.017453292519943295],
AXIS["Geodetic longitude", EAST],
AXIS["Geodetic latitude", NORTH]],
PROJECTION["Mercator_1SP"],
PARAMETER["latitude_of_origin", 0.0],
PARAMETER["central_meridian", 0.0],
PARAMETER["scale_factor", 1.0],
PARAMETER["false_easting", 0.0],
PARAMETER["false_northing", 0.0],
UNIT["m", 1.0],
AXIS["Easting", EAST],
AXIS["Northing", NORTH]]</crs>
<resampling>Nearest</resampling>
<referencePixelX/>
<referencePixelY/>
<easting/>
<northing/>
<orientation/>
<pixelSizeX/>
<pixelSizeY/>
<width/>
<height/>
<tileSizeX/>
<tileSizeY/>
<orthorectify>false</orthorectify>
<elevationModelName/>
<noDataValue>NaN</noDataValue>
<includeTiePointGrids>true</includeTiePointGrids>
<addDeltaBands>false</addDeltaBands>
</parameters>
</node>
<node id="Subset">
<operator>Subset</operator>
<sources>
<sourceProduct refid="Reproject"/>
</sources>
<parameters class="com.bc.ceres.binding.dom.XppDomElement">
<sourceBands/>
<region>0,0,0,0</region>
<referenceBand/>
<geoRegion>POLYGON ((${east} ${north}, ${west} ${north}, ${west} ${south}, ${east} ${south}, ${east} ${north}))</geoRegion>
<subSamplingX>1</subSamplingX>
<subSamplingY>1</subSamplingY>
<fullSwath>false</fullSwath>
<tiePointGridNames/>
<copyMetadata>true</copyMetadata>
</parameters>
</node>
<node id="Write">
<operator>Write</operator>
<sources>
<source>Subset</source>
</sources>
<parameters>
<file>${ofile}</file>
<formatName>NetCDF4-CF</formatName>
</parameters>
</node>
</graph>
------------------------------------------------------
*file.par*
east=-51
west=-63
north=-35
south=-42
referenceComponentPattern=${ORIGINAL_NAME}_chl
dependentComponentPattern=${ORIGINAL_NAME}_sst