Re: imagery data: geometry shift validation

# -*- coding: utf-8 -*-
"""
Created on Fri Feb 25 2015

@author: Sajid Pareeth
@email: [email protected]
"""
import os
#set following variables first
os.environ['PPP_CONFIG_DIR'] = '/usr/local/src/mpop/etc'
os.environ['PYGAC_CONFIG_FILE'] = '/usr/local/src/pygac/etc/pygac.cfg'
# Processing of avhrr l1b bands to tiff starts here...
import numpy as np
import pyresample
import pyproj
import geotiepoints
import argparse
import datetime
from pyresample import *
from pyresample.gradient_search import gradient_search
from pyresample.geometry import SwathDefinition
from pyresample.geometry import AreaDefinition
from mpop.satellites import PolarFactory
from mpop.projector import get_area_def
from mpop.satin.lac_klm_l1b import KLMReader

def main():
    usage = "usage: %prog [arguments]"
    parser = argparse.ArgumentParser(description='Process AVHRR L1B LAC KLM data from NOAA-CLASS')
    parser.add_argument('-n', "--name", dest="name",
                   help='name of the image file (L1B)')
    parser.add_argument('-I', "--instrument", dest="instrument",
                   help='NOAA Instrument code - NL(NOAA16),NM(NOAA17),NN(NOAA18),NP(NOAA19)')
    parser.add_argument('-Y', "--year", dest="year",
                   help='Year of the image')
    parser.add_argument('-d', "--doy", dest="doy",
                   help='Day of the year of the image')
    parser.add_argument('-H', "--hour", dest="hour",
                   help='Hour of the image in 24 hour format')
    parser.add_argument('-M', "--minute", dest="minute",
                   help='Minute of the image')
    parser.add_argument('-i', "--ipdirectory", dest="ipdirectory",
                   help='path to folder where input images are stored')
    parser.add_argument('-o', "--opdirectory", dest="opdirectory",
                   help='path to folder where output are saved')
		   
    args = parser.parse_args()
    mname = args.name
    minst = args.instrument
    myear = args.year
    mdoy = args.doy
    mhour = args.hour
    mminute = args.minute
    mipdir = args.ipdirectory
    mopdir = args.opdirectory
    #Pytroll starts here    
    reader = KLMReader()
    ts = datetime.datetime.strptime('{year} {doy} {hour} {minute}'.format(year=myear,doy=mdoy,hour=mhour,minute=mminute), '%y %j %H %M')
    reader.read("{ipdir}/{name}".format(ipdir=mipdir,name=mname))
    reader.get_lonlat()
    area_swath = SwathDefinition(reader.lons, reader.lats)
    channels = reader.get_calibrated_channels()
    sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles()
    ID=reader.instrument_id
    global_data = PolarFactory.create_scene("noaa", "{ID}".format(ID=ID), "avhrr", ts)
    area_def = get_area_def("euro_laea_AVHRR")
    scene = global_data
    global_data.area = area_swath
    # To save zenith angle as tiff, use the following line...
    channels[:, :, 2] = sat_zen
    scene[0.63] = channels[:, :, 0]
    scene[0.9125] = channels[:, :, 1]
    scene[3.74] = channels[:, :, 2]
    scene[10.8] = channels[:, :, 3]
    scene[12.0] = channels[:, :, 4]
    scene.area = area_swath
    l = scene.project("euro_laea_AVHRR")
    l1 = gradient_search(scene[0.63].data.astype(np.float64),scene[0.63].area.lons,scene[0.63].area.lats,area_def)
    l[0.63] = np.ma.masked_values(l1, 0)
    b1 = l.image(0.63, mode="L")
    b1.save("{opdir}/NSS.LHRR.{inst}.D{year}{doy}.S{hour}{minute}_b1.tif".format(opdir=mopdir,inst=minst,year=myear,doy=mdoy,hour=mhour,minute=mminute), floating_point=True)
    l2 = gradient_search(scene[0.9125].data.astype(np.float64),scene[0.9125].area.lons,scene[0.9125].area.lats,area_def)
    l[0.9125] = np.ma.masked_values(l2, 0)
    b2 = l.image(0.9125, mode="L")
    b2.save("{opdir}/NSS.LHRR.{inst}.D{year}{doy}.S{hour}{minute}_b2.tif".format(opdir=mopdir,inst=minst,year=myear,doy=mdoy,hour=mhour,minute=mminute), floating_point=True)
    lz = gradient_search(scene[3.74].data.astype(np.float64),scene[3.74].area.lons,scene[3.74].area.lats,area_def)
    l[3.74] = np.ma.masked_values(lz, 0)
    bz = l.image(3.74, mode="L")
    bz.save("{opdir}/NSS.LHRR.{inst}.D{year}{doy}.S{hour}{minute}_bz.tif".format(opdir=mopdir,inst=minst,year=myear,doy=mdoy,hour=mhour,minute=mminute), floating_point=True)
    l4 = gradient_search(scene[10.8].data.astype(np.float64),scene[10.8].area.lons,scene[10.8].area.lats,area_def)
    l[10.8] = np.ma.masked_values(l4, 0)
    b4 = l.image(10.8, mode="L")
    b4.save("{opdir}/NSS.LHRR.{inst}.D{year}{doy}.S{hour}{minute}_b4.tif".format(opdir=mopdir,inst=minst,year=myear,doy=mdoy,hour=mhour,minute=mminute), floating_point=True)
    l5 = gradient_search(scene[12.0].data.astype(np.float64),scene[12.0].area.lons,scene[12.0].area.lats,area_def)
    l[12.0] = np.ma.masked_values(l5, 0)
    b5 = l.image(12.0, mode="L")
    b5.save("{opdir}/NSS.LHRR.{inst}.D{year}{doy}.S{hour}{minute}_b5.tif".format(opdir=mopdir,inst=minst,year=myear,doy=mdoy,hour=mhour,minute=mminute), floating_point=True)
   
    ################ENDS HERE############################################

if __name__ == "__main__":
    main()

#!/bin/sh
## Processing NOAA AVHRR LAC L1B images###
## Author: Sajid Pareeth, 2015
##The code below will not run as such, you have to adapt it to the local environment
##The below code should run inside GRASS GIS session

##Exit strategy if you are not inside GRASS GIS
if [ -z "$GISBASE" ] ; then
    echo "You must be in GRASS GIS to run this program." >&2
    exit 1
fi
##Setting the GRASS environments
export GRASS_OVERWRITE=1
export GRASS_MESSAGE_FORMAT=plain  # percent output as 0..1..2..
# setting environment, so that awk works properly in all languages
unset LC_ALL
LC_NUMERIC=C
export LC_NUMERIC

#Setting the study area region
g.region n=2531000 s=2480000 w=4360000 e=4390000 res=$RES -a

#setting the Python environment:
PYVERSION=`python --version 2>&1 | cut -d' ' -f2 | cut -d'.' -f1-2`
MYPYSITES=/usr/local/lib64/python$PYVERSION
# for runtime, attach to existing PYTHONPATH:
export PYTHONPATH=$PYTHONPATH:$MYPYSITES/site-packages
#setting paths required for PYTROLL
PPP_CONFIG_DIR=/home/sajid/mpop_etc
export PPP_CONFIG_DIR=$PPP_CONFIG_DIR
PYGAC_CONFIG_FILE=/home/sajid/pygac_etc/pygac.cfg
export PYGAC_CONFIG_FILE=$PYGAC_CONFIG_FILE

#The loop over all the days from 1986 to 2014 starts here
for yyyy in `seq 1986 2014`; do
    cd ${MYDATA}
    yy=`echo ${yyyy}|cut -c3-4`
    leap=`is_leap_year.sh $yyyy` #small script to check leap year or not
    if [ $leap -eq 1 ]; then
	nd=366
    else
	nd=365
    fi
    for d in `seq 1 $nd`; do
     	cd ${MYDATA}
	dy=`echo $d | awk '{ printf("%03d\n", $1) }'`
	NUM=`ls NSS.LHRR.NP.D${yy}${dy}*|wc -l`
	if [ ${NUM} -eq 0 ]; then
	    echo "No images on ${yyyy}${dy}"
	    continue
	else
	    for y in `ls NSS.LHRR.NP.D${yy}${dy}*`; do
		echo "##########Processing ${y} starts here#########"
		cd ${MYDATA}/"LACdata"${yyyy}"_NOheader"
		# Set the region
		g.region n=$N s=$S w=$W e=$E res=$RES -a
		YEAR=`echo ${y}|cut -c14-15`
		DOY=`echo $y|cut -c16-18`
		HOUR=`echo $y|cut -c21-22`
		MIN=`echo $y|cut -c23-24`
		TIME=`echo $y|cut -c21-24`
		# i=NSS.LHRR.NP.D14324.S1121
		i=`echo $y|cut -c1-24`
		# j=NSS.LHRR.NP.D14324
		j=`echo $y|cut -c1-18`
		if [ ${TIME} -ge ${MINTIME} -a ${TIME} -le ${MAXTIME} ]; then
		    echo "##PYTROLL starts here##"
		    python pytroll.py -i ${y} -o ${OUTDIR}
		    echo "##PYTROLL ends here##"
		    #Geo-correction using OTBcli_homologous points
		    otbcli_HomologousPointsExtraction -in1 input_b1.tif -band1 1 -in2 ref_b1.tif -band2 1 -algorithm sift -mode full -out $OUTB1
		    otbcli_HomologousPointsExtraction -in1 input_b1.tif -band1 1 -in2 ref_b1.tif -band2 1 -algorithm sift -mode full -out $OUTB2
                    otbcli_HomologousPointsExtraction -in1 input_b1.tif -band1 1 -in2 ref_b1.tif -band2 1 -algorithm sift -mode full -out $OUTB3
		    otbcli_HomologousPointsExtraction -in1 input_b1.tif -band1 1 -in2 ref_b1.tif -band2 1 -algorithm sift -mode full -out $OUTB4
                    echo "##Feature matchin (SIFT) using OTB ends here##"
		    ##To make the tie-points from SIFT Grass compatible
		    cat $OUTB1 $OUTB2 $OUTB3 $OUTB4 > ${MYTMPDIR}/${i}_all.txt 
		    ## adding a column with enable/disable
		    awk '{$5="1\t"$5}1' ${MYTMPDIR}/${i}_all.txt > ${MYTMPDIR}/${i}_all_GRASS.txt 
		    # proceed in the GRASS GIS database
		    # importing the input TIFF files
		    r.in.gdal input=${y}_b1.tif output=${i}_b1 memory=${MEMORY}
		    r.in.gdal input=${y}_b2.tif output=${i}_b2 memory=${MEMORY}
		    r.in.gdal input=${y}_bz.tif output=${i}_zenith memory=${MEMORY}
		    r.in.gdal input=${y}_b4.tif output=${i}_b4 memory=${MEMORY}
		    r.in.gdal input=${y}_b5.tif output=${i}_b5 memory=${MEMORY}
		    i.group group=${y} input=${y}_b1,${y}_b2,${y}_b4,${y}_b5,${y}_zenith
		    #Moving the POINTS file to the GRASS group folder
		    mv ${MYTMPDIR}/${i}_all_GRASS.txt ${GRASSLOC}/group/${i}/POINTS
		    ##GCP filtering and geo-rectification
		    PTCNT=`wc -l ${GRASSLOC}/group/${i}/POINTS|cut -d' ' -f1`
		    unset use 
		    if [ ${PTCNT} -le 3 ]; then
			g.remove group name=${i} -f
			echo "Not enough gcps to filter (< 3), hence ${y} is not usable"
			continue
		    else
			eval `m.gcp.filter group=${i} order=1 threshold=500 -b`
			USE=${use}
			i.target group=${i} -c
		    fi
		    
		    if [ ${USE} -lt 20 ]; then
			#Exit, delete all the above images and go to the next one
			#cleanup "$i"
			g.remove group name=${i} -f
			echo "${y} is not usable due to lack of enough homologous points - ${USE}, hence avoiding"
			continue
		    elif [ ${USE} -ge 20 ] && [ ${USE} -le 300 ]; then
			m.gcp.filter group=${i} order=1 threshold=500 -b -u
			i.rectify -a group=${i} extension=_rectified order=1 method=nearest --o
		    elif [ ${USE} -gt 300 ]; then
			m.gcp.filter group=${i} order=2 threshold=500 -b -u
			i.rectify -a group=${i} extension=_rectified order=2 method=nearest --o
		    fi
		    unset use
		    r.colors map=${i}_b1_rectified rules=${CLRDIR}/avhrrb1.clr
		    r.colors map=${i}_b2_rectified rules=${CLRDIR}/avhrrb2.clr
		    r.colors map=${i}_b4_rectified color=kelvin
		    r.colors map=${i}_b5_rectified color=kelvin
		    #Cloud mask based on SPARC 
		    echo "SPARC Cloud detection starts here"
		    if [ ${MIN} -lt 30 ]; then
			ECMWFHR=${HOUR}
		    elif [ ${MIN} -gt 30 ] && [ ${HOUR} -eq 23 ]; then
			ECMWFHR=00
		    else
			ECMWFHR=`echo $((${HOUR} + 1))`
		    fi
		    
		    D=`echo $((${d} - 1))`
		    DATE=`date -d "${D} days ${yyyy}-01-01" +"%d%m%Y"`
		    M=00
		    echo "${yyyy}${DOY}_${ECMWFHR}${M}_${DATE}_ecmwf@sp_ecmwf is used for cloud removal"  
		    #2006252_1200_09092006_ecmwf - sample ecmwf data
		    #below the offset and scale factors from Trischenko et.al 2006
		    r.mapcalc "${i}_Tindex = 1.0 * (${i}_b4_rectified - ${yyyy}${DOY}_${ECMWFHR}${M}_${DATE}_ecmwf_hants@sp_ecmwf + 6.0) * -0.42"
		    r.mapcalc "${i}_Bindex_land = 1.0 * (${i}_b1_rectified - 0.30) * 67"
		    r.mapcalc "${i}_Bindex_water = 1.0 * (${i}_b2_rectified - 0.24) * 75"
		    r.mapcalc "${i}_Cindex = 1.0 * (${i}_b4_rectified - ${i}_b5_rectified - 1.5) * 4.0"
		    #Please try with Bindex condition to B < -12 instead of B > -12	in case of problematic index values
		    #B_test_land     = if(${i}_Bindex_land < -12 && isnull(watermask), null(), T_test), 
		    #B_test_water    = if(${i}_Bindex_water < -12 && watermask == 1, null(), B_test_land),
		    r.mapcalc << EOF
		    ${i}_b4_masked = eval( \\
		    T_test          = if(${i}_Tindex > 8, null(), ${i}_b4_rectified), \\
		    C_test          = if(${i}_Cindex > 25, null(), T_test), \\
		    viewangle_test  = if(${i}_zenith_rectified > 45, null(), C_test), \\
		    viewangle_test)
EOF
		    r.mapcalc << EOF
		    ${i}_b5_masked = eval( \\
		    T_test          = if(${i}_Tindex > 8, null(), ${i}_b5_rectified), \\
		    C_test          = if(${i}_Cindex > 25, null(), T_test), \\
		    viewangle_test  = if(${i}_zenith_rectified > 45, null(), C_test), \\
		    viewangle_test)
EOF
		    r.mapcalc "${i}_zenith_masked = if(isnull(${i}_b4_masked), null(), ${i}_zenith_rectified)"
		    echo "SPARC Cloud detection ends here"
		else
		    echo "${y} is outside the time frame, hence avoiding";
		    continue
		fi
	    done
	fi
    done
    echo "#################Processing ${yyyy} finishes here#############"
done