From 05ed08ca28bdb7fd0baef2ea0dbccb5e057c9dee Mon Sep 17 00:00:00 2001 From: Dale Roberts Date: Mon, 6 Jun 2016 11:46:55 +1000 Subject: [PATCH 1/5] Vector version of surface_reflectance.f90 (and called functions) discussed in meeting on Friday 3/6/16. Also included are changes in the makefiles necessary for an Intel compiler build, and updated nbar.cfg with optimal y_tile size. --- gaip/reflectance_mod.f90 | 99 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 99 insertions(+) create mode 100644 gaip/reflectance_mod.f90 diff --git a/gaip/reflectance_mod.f90 b/gaip/reflectance_mod.f90 new file mode 100644 index 0000000..d04bf5a --- /dev/null +++ b/gaip/reflectance_mod.f90 @@ -0,0 +1,99 @@ +module reflectance_mod + + real, parameter :: pi = 3.14159265358979 + real, parameter :: pib = 0.017453292519943 ! pi / 180 + +contains + + function RL_brdf (solar,view,ra,hb,br,brdf0,brdf1,brdf2, veclen) +! calculate the normalised BRDF shape function + integer :: veclen + real solar(veclen), view(veclen), ra(veclen), hb, br + real brdf0,brdf1,brdf2 + real RL_brdf(veclen) + + real rs_thick(veclen),li_sparse(veclen),secsolar(veclen),secvia(veclen) + real cossolar(veclen),cosvia(veclen),cosra(veclen),sinsolar(veclen),sinvia(veclen),sinra(veclen) + real cosxi(veclen),xi(veclen),tansolar(veclen),tanvia(veclen),theta_new_v(veclen),theta_new_s(veclen) + real d_li2(veclen),x_li(veclen),cosl(veclen),l_li(veclen),o_li(veclen) + + real tan_theta_new_v(veclen), cos_theta_new_v(veclen), sin_theta_new_v(veclen) + real tan_theta_new_s(veclen), cos_theta_new_s(veclen), sin_theta_new_s(veclen) + + RL_brdf=0.0 + +! pi = 4 * atan(1.0) +! calculate Ross-thick kernel + cossolar = cos(solar) + cosvia = cos(view) + cosra = cos(ra) + sinsolar = sin(solar) + sinvia = sin(view) + sinra = sin(ra) + cosxi = min(cossolar*cosvia+sinsolar*sinvia*cosra,1.0) +! if (cosxi .ge.1) then +! cosxi = 1 +! endif + xi = acos(cosxi) + rs_thick = ((pi/2.0-xi)*cosxi+sin(xi)) / (cossolar+cosvia) - pi/4.0 + +! alculate Li-sparse + tansolar = sinsolar/cossolar + tanvia = sinvia/cosvia +! theta_new_v = atan(br*tanvia) +! theta_new_s = atan(br*tansolar) + tan_theta_new_v = br*tanvia + ! Trig identities + cos_theta_new_v = 1.0 / sqrt(1.0 + tan_theta_new_v * tan_theta_new_v) + sin_theta_new_v = tan_theta_new_v * cos_theta_new_v + + tan_theta_new_s = br*tansolar + ! Trig identities + cos_theta_new_s = 1.0 / sqrt(1.0 + tan_theta_new_s * tan_theta_new_s) + sin_theta_new_s = tan_theta_new_s * cos_theta_new_s + +! cosxi = cos(theta_new_s)*cos(theta_new_v)+sin(theta_new_s)*sin(theta_new_v)*cosra + cosxi = min( cos_theta_new_s*cos_theta_new_v + sin_theta_new_s*sin_theta_new_v*cosra, 1.0 ) +! if (cosxi.ge.1) then +! cosxi = 1 +! endif + secsolar = 1.0/cos_theta_new_s + secvia = 1.0/cos_theta_new_v +! d_li2 = abs(tan(theta_new_s)**2+tan(theta_new_v)**2 - 2.0*tan(theta_new_s)*tan(theta_new_v)*cosra) + d_li2 = abs(tan_theta_new_s**2+tan_theta_new_v**2 - 2.0*tan_theta_new_s*tan_theta_new_v*cosra) + x_li = tan_theta_new_s*tan_theta_new_v*sinra + cosl = hb*sqrt(d_li2+x_li**2)/(secsolar+secvia) +! if (cosl .ge. 1.0) then +! o_li=0.0 +! else + l_li=acos(cosl) + o_li=(l_li-sin(l_li)*cosl)*(secsolar+secvia)/pi + where( cosl .ge. 1.0 ) o_li = 0.0 +! endif + li_sparse=o_li-(secsolar+secvia)+0.5*(1.0+cosxi) * secsolar*secvia + RL_brdf=brdf0+brdf1*rs_thick+brdf2*li_sparse + return + end function RL_brdf + + function black_sky(brdf0,brdf1,brdf2,theta,veclen) + integer veclen + real brdf0,brdf1,brdf2 + real theta(veclen) + real black_sky(veclen) + +! theta should be in radians + black_sky = brdf0 + & + brdf1*(-0.007574-0.070987*theta**2 + 0.307588*theta**3) + & + brdf2*(-1.284909-0.166314*theta**2 + 0.041840*theta**3) + return + end function black_sky + + + function white_sky(brdf0, brdf1, brdf2) + real brdf0, brdf1, brdf2 + real white_sky + white_sky = brdf0 + brdf1*0.189184 - brdf2*1.377622 + return + end function white_sky + +end module reflectance_mod From 0b298ba4701dad0d2111b08af178badca4ed733b Mon Sep 17 00:00:00 2001 From: Dale Roberts Date: Mon, 6 Jun 2016 12:13:05 +1000 Subject: [PATCH 2/5] Updated working copy to latest revision, merged changes --- .gitignore | 1 + gaip/Makefile | 4 +- gaip/black_sky.f90 | 6 +- gaip/brdf_shape.f90 | 82 ++++-- gaip/setup.py | 4 +- gaip/surface_reflectance.f90 | 537 +++++++++++++++++++++++------------ src/Makefile | 6 +- workflow/nbar.cfg | 2 +- 8 files changed, 414 insertions(+), 228 deletions(-) diff --git a/.gitignore b/.gitignore index 6df902f..b72b97c 100644 --- a/.gitignore +++ b/.gitignore @@ -35,3 +35,4 @@ nosetests.xml .idea *.iml *.ipr +*~ \ No newline at end of file diff --git a/gaip/Makefile b/gaip/Makefile index 3be082d..be5723d 100644 --- a/gaip/Makefile +++ b/gaip/Makefile @@ -1,4 +1,4 @@ -FC := gfortran +FC := ifort F2PY:= f2py all: _cast_shadow_mask.so _slope_self_shadow.so _surface_reflectance.so sys_variables.o sys_variables.mod _sat_sol_angles.so _satellite_model.so _track_time_info.so unmiximage.pyf unmiximage.so _interpolation.so @@ -6,7 +6,7 @@ all: _cast_shadow_mask.so _slope_self_shadow.so _surface_reflectance.so sys_vari _cast_shadow_mask.so _slope_self_shadow.so _surface_reflectance.so sys_variables.o sys_variables.mod _sat_sol_angles.so _satellite_model.so _track_time_info.so: $$(module load gcc/5.2.0) $(FC) -c sys_variables.f90 - python setup.py build_ext -i + python setup.py config_fc --opt "-O3 -g" --arch "-xHost" build_ext -i unmiximage.pyf unmiximage.so: $(F2PY) -h unmiximage.pyf -m unmiximage unmiximage.f90 diff --git a/gaip/black_sky.f90 b/gaip/black_sky.f90 index 4168623..40f3bbe 100644 --- a/gaip/black_sky.f90 +++ b/gaip/black_sky.f90 @@ -1,6 +1,8 @@ -real function black_sky(brdf0,brdf1,brdf2,theta) +function black_sky(brdf0,brdf1,brdf2,theta,veclen) + integer veclen real brdf0,brdf1,brdf2 - real theta + real theta(veclen) + real black_sky(veclen) ! theta should be in radians black_sky = brdf0 + & diff --git a/gaip/brdf_shape.f90 b/gaip/brdf_shape.f90 index 76aeba0..a36c33a 100644 --- a/gaip/brdf_shape.f90 +++ b/gaip/brdf_shape.f90 @@ -1,16 +1,23 @@ -real function RL_brdf (solar,view,ra,hb,br,brdf0,brdf1,brdf2) +function RL_brdf (solar,view,ra,hb,br,brdf0,brdf1,brdf2,veclen) ! calculate the normalised BRDF shape function - real pi - real solar, view, ra, hb, br + integer veclen + real solar(veclen), view(veclen), ra(veclen), hb, br real brdf0,brdf1,brdf2 - real rs_thick,li_sparse,secsolar,secvia - real cossolar,cosvia,cosra,sinsolar,sinvia,sinra - real cosxi,xi,tansolar,tanvia,theta_new_v,theta_new_s - real d_li2,x_li,cosl,l_li,o_li + real RL_brdf(veclen) + + real rs_thick(veclen),li_sparse(veclen),secsolar(veclen),secvia(veclen) + real cossolar(veclen),cosvia(veclen),cosra(veclen),sinsolar(veclen),sinvia(veclen),sinra(veclen) + real cosxi(veclen),xi(veclen),tansolar(veclen),tanvia(veclen),theta_new_v(veclen),theta_new_s(veclen) + real d_li2(veclen),x_li(veclen),cosl(veclen),l_li(veclen),o_li(veclen) + + real tan_theta_new_v(veclen), cos_theta_new_v(veclen), sin_theta_new_v(veclen) + real tan_theta_new_s(veclen), cos_theta_new_s(veclen), sin_theta_new_s(veclen) + + real, parameter :: pi = 3.14159265358979 RL_brdf=0.0 - pi = 4 * atan(1.0) +! pi = 4 * atan(1.0) ! calculate Ross-thick kernel cossolar = cos(solar) cosvia = cos(view) @@ -18,34 +25,47 @@ real function RL_brdf (solar,view,ra,hb,br,brdf0,brdf1,brdf2) sinsolar = sin(solar) sinvia = sin(view) sinra = sin(ra) - cosxi = cossolar*cosvia+sinsolar*sinvia*cosra - if (cosxi .ge.1) then - cosxi = 1 - endif + cosxi = min(cossolar*cosvia+sinsolar*sinvia*cosra,1.0) +! if (cosxi .ge.1) then +! cosxi = 1 +! endif xi = acos(cosxi) - rs_thick = ((pi/2.0-xi)*cos(xi)+sin(xi)) / (cossolar+cosvia) - pi/4.0 + rs_thick = ((pi/2.0-xi)*cosxi+sin(xi)) / (cossolar+cosvia) - pi/4.0 ! alculate Li-sparse tansolar = sinsolar/cossolar tanvia = sinvia/cosvia - theta_new_v = atan(br*tanvia) - theta_new_s = atan(br*tansolar) - cosxi = cos(theta_new_s)*cos(theta_new_v)+sin(theta_new_s)*sin(theta_new_v)*cosra - if (cosxi.ge.1) then - cosxi = 1 - endif - secsolar = 1.0/cos(theta_new_s) - secvia = 1.0/cos(theta_new_v) - d_li2 = abs(tan(theta_new_s)**2+tan(theta_new_v)**2 - 2.0*tan(theta_new_s)*tan(theta_new_v)*cosra) - x_li = tan(theta_new_s)*tan(theta_new_v)*sinra +! theta_new_v = atan(br*tanvia) +! theta_new_s = atan(br*tansolar) + tan_theta_new_v = br*tanvia + ! Trig identities + cos_theta_new_v = 1.0 / sqrt(1 + tan_theta_new_v * tan_theta_new_v) + sin_theta_new_v = tan_theta_new_v * cos_theta_new_v + + tan_theta_new_s = br*tansolar + ! Trig identities + cos_theta_new_s = 1.0 / sqrt(1 + tan_theta_new_s * tan_theta_new_s) + sin_theta_new_s = tan_theta_new_s * cos_theta_new_s + +! cosxi = cos(theta_new_s)*cos(theta_new_v)+sin(theta_new_s)*sin(theta_new_v)*cosra + cosxi = min( cos_theta_new_s*cos_theta_new_v + sin_theta_new_s*sin_theta_new_s*cosra, 1.0 ) +! if (cosxi.ge.1) then +! cosxi = 1 +! endif + secsolar = 1.0/cos_theta_new_s + secvia = 1.0/cos_theta_new_v +! d_li2 = abs(tan(theta_new_s)**2+tan(theta_new_v)**2 - 2.0*tan(theta_new_s)*tan(theta_new_v)*cosra) + d_li2 = abs(tan_theta_new_s**2+tan_theta_new_v**2 - 2.0*tan_theta_new_s*tan_theta_new_v*cosra) + x_li = tan_theta_new_s*tan_theta_new_v*sinra cosl = hb*sqrt(d_li2+x_li**2)/(secsolar+secvia) - if (cosl .ge. 1.0) then - o_li=0.0 - else - l_li=acos(cosl) - o_li=(l_li-sin(l_li)*cos(l_li))*(secsolar+secvia)/pi - endif - li_sparse=o_li-(secsolar+secvia)+0.5*(1.0+cosxi) * secsolar*secvia - RL_brdf=brdf0+brdf1*rs_thick+brdf2*li_sparse +! if (cosl .ge. 1.0) then +! o_li=0.0 +! else + l_li=acos(cosl) + o_li=(l_li-sin(l_li)*cosl)*(secsolar+secvia)/pi + where( cosl .ge. 1.0 ) o_li = 0.0 +! endif + li_sparse=o_li-(secsolar+secvia)+0.5*(1.0+cosxi) * secsolar*secvia + RL_brdf=brdf0+brdf1*rs_thick+brdf2*li_sparse return end diff --git a/gaip/setup.py b/gaip/setup.py index 19fd653..a9d2682 100644 --- a/gaip/setup.py +++ b/gaip/setup.py @@ -25,9 +25,7 @@ 'geo2metres_pixel_size.f90']), Extension(name='_surface_reflectance', sources=['surface_reflectance.f90', - 'white_sky.f90', - 'black_sky.f90', - 'brdf_shape.f90']), + 'reflectance_mod.f90']), Extension(name='_satellite_model', sources=['geo2metres_pixel_size.f90', 'satellite_model.f90']), diff --git a/gaip/surface_reflectance.f90 b/gaip/surface_reflectance.f90 index 4b14793..2b6e0bc 100644 --- a/gaip/surface_reflectance.f90 +++ b/gaip/surface_reflectance.f90 @@ -36,7 +36,9 @@ SUBROUTINE reflectance( & ! Calculates lambertian, brdf corrected and terrain corrected surface ! reflectance. + use reflectance_mod + implicit none ! input parameters integer nrow, ncol ! we should be passing in transposed arrays cols = rows and vice versa real*4 rori ! threshold for terrain correction @@ -89,29 +91,61 @@ SUBROUTINE reflectance( & !f2py intent(inout) iref_lm, iref_brdf, iref_terrain ! internal parameters + real, parameter :: hb = 2.0 + real, parameter :: br = 1.0 + integer, parameter :: veclen = 16 + + !!! Tiled input quantities + real*4 tile_dn_1(veclen) ! raw image (in real4 as replacement for dn array) + integer*2 tile_mask_self(veclen) ! mask + integer*2 tile_mask_castsun(veclen) ! self shadow mask + integer*2 tile_mask_castview(veclen) ! cast shadow mask +! real*4 tile_solar_angle(veclen) ! solar zenith angle +! real*4 tile_sazi_angle(veclen) ! solar azimuth angle +! real*4 tile_view_angle(veclen) ! view angle (for flat surface) + real*4 tile_rela_angle(veclen) ! relative azimuth angle (for flat surface) +! real*4 tile_slope_angle(veclen) ! slop angle +! real*4 tile_aspect_angle(veclen) ! aspect angle +! real*4 tile_it_angle(veclen) ! incident angle (for inclined surface) +! real*4 tile_et_angle(veclen) ! exiting angle (for inclined surface) + real*4 tile_rela_slope(veclen) ! relative angle (for inclined surface) + real*4 tile_a_mod(veclen) ! modtran output (a) + real*4 tile_b_mod(veclen) ! modtran output (b) + real*4 tile_s_mod(veclen) ! modtran output (s) + real*4 tile_fv(veclen) ! modtran output (fs) + real*4 tile_fs(veclen) ! modtran output (fv) + real*4 tile_ts(veclen) ! modtran output (ts) + real*4 tile_edir_h(veclen) ! modtran output (direct irradiance) + real*4 tile_edif_h(veclen) ! modtran output (diffuse irradiance) + integer*2 tile_iref_lm(veclen) ! atmospheric corrected lambertial reflectance + integer*2 tile_iref_brdf(veclen) ! atmospheric and brdf corrected reflectance + integer*2 tile_iref_terrain(veclen) ! atmospheric and brdf and terrain correc + + real*4 tile_ref_lm(veclen) + real*4 tile_ref_brdf(veclen) + real*4 tile_ref_terrain(veclen) +!!! Tile sizing + integer :: istart, iend, tilelen + integer i, j - real pi, pib - real ann_f, aa_viewf, aa_solarf, aa_white - real ann_s, aa_views, aa_solars - real lt, norm_1, norm_2 - real solar, sazi, view, slope, aspect, ra_lm, ra_sl, it, et - - double precision a_eqf, b_eqf, c_eqf, ref_barf, aa_flat - double precision a_eqs, b_eqs, c_eqs, ref_bars, aa_slope - double precision ref_brdfrealf, ref_brdfreals - double precision bb_angle, cc_angle, tttt - real hb, br, fnn, it_brdf, it_bk, et_brdf, et_bk - real angle_th - real vt, vd,angle, edir_t, eadj, edif_t, rdir, rdif, rtotal - real RL_brdf, black_sky, white_sky - real cosslope, rth - external RL_brdf, black_sky, white_sky + real ann_f(veclen), aa_viewf(veclen), aa_solarf(veclen), aa_white + real ann_s(veclen), aa_views(veclen), aa_solars(veclen) + real lt(veclen), norm_1, norm_2 + real solar(veclen), sazi(veclen), view(veclen), slope(veclen), aspect(veclen), ra_lm(veclen), ra_sl(veclen), it(veclen), et(veclen) + + double precision a_eqf(veclen), b_eqf(veclen), c_eqf, ref_barf(veclen), aa_flat(veclen) + double precision a_eqs(veclen), b_eqs(veclen), c_eqs, ref_bars(veclen), aa_slope(veclen) + double precision ref_brdfrealf(veclen), ref_brdfreals(veclen) + double precision bb_angle(veclen), cc_angle(veclen), tttt(veclen) + real fnn, fnn_a(1), it_brdf(veclen), it_bk(veclen), et_brdf(veclen), et_bk(veclen) + real angle_th(veclen) + real vt(veclen), vd(veclen),angle(veclen), edir_t(veclen), eadj(veclen), edif_t(veclen), rdir(veclen), rdif(veclen), rtotal(veclen), rtotal_tmp(veclen), rdir_tmp(veclen) + real cossolar(veclen),cosit(veclen),rth(veclen) + + logical :: tilemask(veclen) ! li-sparse parameters - hb = 2.0 - br = 1.0 - pi = 4 * atan(1.0) - pib = pi / 180.0 + norm_1 = brdf1 / brdf0 norm_2 = brdf2 / brdf0 @@ -119,228 +153,359 @@ SUBROUTINE reflectance( & ! calculate white sky albedo aa_white = white_sky(1.0, norm_1, norm_2) ! calcualte BRDF at 45 solar angle and 0 view angle - fnn = RL_brdf(45 * pib, 0.0, 0.0, hb, br, 1.0, norm_1, norm_2) -! print*,fnn + fnn_a = RL_brdf( (/ 45 * pib /), (/ 0.0 /), (/ 0.0 /), hb, br, 1.0, norm_1, norm_2,1) + fnn = fnn_a(1) ! Now loop over the cols of the images do j=1,ncol !---------------------------------------------------------------------- ! convert byte to integer for the raw image - do i=1,nrow - if (dn_1(i, j) .lt. 0) then - dn(i) = dn_1(i, j) + 65536 - else - dn(i) = dn_1(i, j) - endif - enddo +! do i=1,nrow +! if (dn_1(i, j) .lt. 0) then +! dn(i) = dn_1(i, j) + 65536 +! else +! dn(i) = dn_1(i, j) +! endif +! enddo + istart = 0 + iend = 0 + do while ( iend < nrow ) + istart = iend + 1 + iend = iend + veclen + if ( iend > nrow ) iend = nrow + tilelen = iend - istart + 1 + + where ( dn_1(istart:iend,j) .lt. 0 ) + tile_dn_1(1:tilelen) = real(dn_1(istart:iend,j) + 65536,4) + elsewhere + tile_dn_1(1:tilelen) = real(dn_1(istart:iend,j),4) + end where + + + if( ANY( a_mod(istart:iend,j) .gt. 0 .AND. tile_dn_1(1:tilelen) .gt. 0.0 ) ) then + + tile_a_mod(1:tilelen) = a_mod(istart:iend,j) + tile_b_mod(1:tilelen) = b_mod(istart:iend,j) + tile_s_mod(1:tilelen) = s_mod(istart:iend,j) + + tile_fs(1:tilelen) = fs(istart:iend,j) + tile_fv(1:tilelen) = fv(istart:iend,j) + tile_ts(1:tilelen) = ts(istart:iend,j) + + tile_edir_h(1:tilelen) = edir_h(istart:iend,j) + tile_edif_h(1:tilelen) = edif_h(istart:iend,j) +!!! Mask tiles + tile_mask_self(1:tilelen) = mask_self(istart:iend,j) + tile_mask_castsun(1:tilelen) = mask_castsun(istart:iend,j) + tile_mask_castview(1:tilelen) = mask_castview(istart:iend,j) + ! now loop over the rows of the images - do i=1,nrow +! do i=1,nrow ! if valid masks and valid digital number then do the calcs - if (a_mod(i, j) .ge. 0 .and. dn(i) .gt. 0) then - if (rela_angle(i, j) .gt. 180) rela_angle(i, j) = rela_angle(i, j) - 360 - if (rela_slope(i, j) .gt. 180) rela_slope(i, j) = rela_slope(i, j) - 360 - +!!! Mask after calculations +! if (a_mod(i, j) .ge. 0 .and. dn(i) .gt. 0) then + +! if (rela_angle(i, j) .gt. 180) rela_angle(i, j) = rela_angle(i, j) - 360 +! if (rela_slope(i, j) .gt. 180) rela_slope(i, j) = rela_slope(i, j) - 360 + where ( rela_angle(istart:iend,j) .gt. 180 ) + tile_rela_angle(1:tilelen) = rela_angle(istart:iend,j) - 360 + elsewhere + tile_rela_angle(1:tilelen) = rela_angle(istart:iend,j) + end where + + where ( rela_slope(istart:iend,j) .gt. 180 ) + tile_rela_slope(1:tilelen) = rela_slope(istart:iend,j) - 360 + elsewhere + tile_rela_slope(1:tilelen) = rela_slope(istart:iend,j) + end where + ! convert angle to radians - solar = solar_angle(i, j) * pib - sazi = sazi_angle(i, j) * pib - view = view_angle(i, j) * pib - slope = slope_angle(i, j) * pib - aspect = aspect_angle(i, j) * pib - ra_lm = rela_angle(i, j) * pib - ra_sl = rela_slope(i, j) * pib - it = it_angle(i, j) * pib - et = et_angle(i, j) * pib - - if (it_angle(i, j) .ge. 70.0) then - it_brdf= 70.0 * pib - else - it_brdf = it_angle(i, j) * pib - endif - - if (it_angle(i, j) .ge. 80.0) then - it_bk = 80.0 * pib - else - it_bk = it_angle(i, j) * pib - endif - - if (et_angle(i, j) .ge. 60.0) then - et_brdf = 60.0 * pib - else - et_brdf = et_angle(i, j) * pib - endif - - if (et_angle(i, j) .ge. 80.0) then - et_bk = 80.0 * pib - else - et_bk = et_angle(i, j) * pib - endif +! solar = solar_angle(i, j) * pib + solar(1:tilelen) = solar_angle(istart:iend,j) + solar = solar * pib + ! Trig is expensive + cossolar = cos(solar) +! sazi = sazi_angle(i, j) * pib + sazi(1:tilelen) = sazi_angle(istart:iend,j) + sazi = sazi * pib + ! view = view_angle(i, j) * pib + view(1:tilelen) = view_angle(istart:iend,j) + view = view * pib + ! slope = slope_angle(i, j) * pib + slope(1:tilelen) = slope_angle(istart:iend,j) + slope = slope * pib + ! aspect = aspect_angle(i, j) * pib + aspect(1:tilelen) = aspect_angle(istart:iend,j) + aspect = aspect * pib + ! ra_lm = rela_angle(i, j) * pib + ra_lm(1:tilelen) = rela_angle(istart:iend,j) + ra_lm = ra_lm * pib + ! ra_sl = rela_slope(i, j) * pib + ra_sl(1:tilelen) = rela_slope(istart:iend,j) + ra_sl = ra_sl * pib + ! it = it_angle(i, j) * pib + it(1:tilelen) = it_angle(istart:iend,j) + ! et = et_angle(i, j) * pib + et(1:tilelen) = et_angle(istart:iend,j) + +! if (it_angle(i, j) .ge. 70.0) then +! it_brdf= 70.0 * pib +! else +! it_brdf = it_angle(i, j) * pib +! endif + it_brdf = min(70.0,it) * pib + +! if (it_angle(i, j) .ge. 80.0) then +! it_bk = 80.0 * pib +! else +! it_bk = it_angle(i, j) * pib +! endif + it_bk = min(80.0, it) * pib + + it = it * pib + ! Trig is expensive + cosit = cos(it) + + ! if (et_angle(i, j) .ge. 60.0) then + ! et_brdf = 60.0 * pib + ! else + ! et_brdf = et_angle(i, j) * pib + ! endif + et_brdf = min(60.0,et) * pib + +! if (et_angle(i, j) .ge. 80.0) then +! et_bk = 80.0 * pib +! else +! et_bk = et_angle(i, j) * pib +! endif + et_bk = min(80.0,et) * pib + + et = et * pib !------------------------------------------------ ! for flat surface ! calculate radiance at top atmosphere - lt = bias + dn(i) * slope_ca +! lt = bias + dn(i) * slope_ca + lt = bias + tile_dn_1 * slope_ca ! calcualte lambetian reflectance with bilinear average - ref_lm(i) = (lt - b_mod(i, j)) / (a_mod(i, j) + s_mod(i, j) * & - (lt - b_mod(i, j))) - iref_lm(i, j) = ref_lm(i) * 10000 + 0.5 +! ref_lm(i) = (lt - b_mod(i, j)) / (a_mod(i, j) + s_mod(i, j) * & +! (lt - b_mod(i, j))) +! iref_lm(i, j) = ref_lm(i) * 10000 + 0.5 + tile_ref_lm = (lt - tile_b_mod) / ( tile_a_mod + tile_s_mod * & + ( lt - tile_b_mod )) + tile_iref_lm = nint(tile_ref_lm * 10000.0,2) ! set as small number if atmospheric corrected reflectance ! below 0.001 - if (ref_lm(i).lt. 0.001) then - ref_lm(i) = 0.001 - ref_brdf(i) = 0.001 - iref_lm(i, j) = 10 - iref_brdf(i, j) = 10 - else +!!! Set these at the end +! if (ref_lm(i).lt. 0.001) then +! ref_lm(i) = 0.001 +! ref_brdf(i) = 0.001 +! iref_lm(i, j) = 10 +! iref_brdf(i, j) = 10 + ! else ! calculate normalized BRDF shape function - ann_f = RL_brdf(solar, view, ra_lm, hb, br, 1.0, norm_1, norm_2) + if( all(tile_ref_lm(1:tilelen) .lt.0.001) ) then + tile_ref_lm = 0.001 + tile_ref_brdf = 0.001 + tile_iref_lm = 10 + tile_iref_brdf = 10 + else + ann_f = RL_brdf(solar, view, ra_lm, hb, br, 1.0, norm_1, norm_2,veclen) ! calculate black sky albedo for sloar angle - aa_solarf = black_sky(1.0, norm_1, norm_2, solar) + aa_solarf = black_sky(1.0, norm_1, norm_2, solar,veclen) ! calculate black sky albedo for view angle - aa_viewf = black_sky(1.0, norm_1, norm_2, view) + aa_viewf = black_sky(1.0, norm_1, norm_2, view,veclen) ! - aa_flat = (fv(i, j) * (fs(i, j) * ann_f + (1.0 - fs(i, j)) * & - aa_viewf) + (1.0 - fv(i, j)) * (fs(i, j) * & - aa_solarf + (1.0 - fs(i, j)) * aa_white)) / aa_white - - a_eqf = (1 - aa_flat) * s_mod(i, j) * (1 - s_mod(i, j) * & - ref_lm(i)) - b_eqf = aa_flat + ref_lm(i) * (1 - aa_flat) * s_mod(i, j) - c_eqf = -ref_lm(i) - - if (abs(a_eqf) .lt. 0.0000001) then - ref_barf = -c_eqf / b_eqf - else - ref_barf = (-b_eqf + & - sqrt(b_eqf**2 - 4 * a_eqf * c_eqf)) / & - (2 * a_eqf) - endif - ref_brdfrealf = ann_f * ref_barf / aa_white - ref_brdf(i) = ref_barf * fnn / aa_white - iref_brdf(i, j) = ref_brdf(i) * 10000 + 0.5 + aa_flat = (tile_fv * (tile_fs * ann_f + (1.0 - tile_fs) * & + aa_viewf) + (1.0 - tile_fv) * (tile_fs * & + aa_solarf + (1.0 - tile_fs) * aa_white)) / aa_white + + a_eqf = (1 - aa_flat) * tile_s_mod * (1 - tile_s_mod * & + tile_ref_lm) + b_eqf = aa_flat + tile_ref_lm * (1 - aa_flat) * tile_s_mod +! c_eqf = -ref_lm(i) + +! if (abs(a_eqf) .lt. 0.0000001) then +! ref_barf = tile_ref_lm / b_eqf +! else + ref_barf = (-b_eqf + & + sqrt(b_eqf**2 + 4 * a_eqf * tile_ref_lm)) / & + (2 * a_eqf) + where( abs(a_eqf) .lt. 0.0000001 ) ref_barf = tile_ref_lm / b_eqf +! endif + ref_brdfrealf = ann_f * ref_barf / aa_white + tile_ref_brdf = min(ref_barf * fnn / aa_white, 1.0) + tile_iref_brdf = nint(tile_ref_brdf * 10000.0,2) +! tile_iref_brdf = int(tile_ref_brdf * 10000.0 + 0.5,2) ! this is to ensure that the brdf correction -! is the same as (or as close as possible to) the original NBAR version - if (ref_brdf(i) .ge. 1) then - ref_brdf(i) = 1.0 - iref_brdf(i, j) = ref_brdf(i) * 10000 - endif - - endif + ! is the same as (or as close as possible to) the original NBAR version +!!! Condition from the top of this section. + where( tile_ref_lm .lt.0.001 ) + tile_ref_lm = 0.001 + tile_ref_brdf = 0.001 + tile_iref_lm = 10 + tile_iref_brdf = 10 + end where + end if + +! if (ref_brdf(i) .ge. 1) then +! ref_brdf(i) = 1.0 +! iref_brdf(i, j) = ref_brdf(i) * 10000 +! endif + +! endif !------------------------------------------------------------------- ! conduct terrain correction - if ((mask_self(i, j) .gt. 0) .and. & - (mask_castsun(i, j).gt. 0) .and. & - (mask_castview(i, j) .gt. 0) .and. (it_angle(i, j) .lt. 90.0) & - .and. (et_angle(i, j) .lt. 90.0)) then +!!! Tile-wide mask at the start, per-element mask at the end + tilemask(1:tilelen) = ((tile_mask_self(1:tilelen) .gt. 0) .and. & + (tile_mask_castsun(1:tilelen) .gt. 0) .and. & + (tile_mask_castview(1:tilelen) .gt. 0) .and. & + (it(1:tilelen) .lt. 85.0*pib ) .and. & + (et(1:tilelen) .lt. 85.0*pib ) ) + if ( ANY( tilemask(1:tilelen) ) ) then !---------------------------------------------------------- - cosslope = cos(slope) ! calculate vd and vt - vd = 0.5 * (1.0 + cosslope) + vd = 0.5 * (1.0 + cos(slope)) vt = 1.0 - vd !--------------------------------------------------------- ! calculate direct irradiance ! Note the account taken of threshold - - edir_t = edir_h(i, j) * cos(it) / cos(solar) + edir_t = tile_edir_h * cosit / cossolar ! calcualte adjacent irradiance for anisotropical surface ! see Iqbal, 1983 "an introduction to solar ! radiation" - eadj = (edir_h(i, j) + edif_h(i, j)) * vt * ref_adj * & - (1.0 + sin(solar / 2.0)**2) * abs(cos(aspect - sazi)) + eadj = (tile_edir_h + tile_edif_h) * vt * ref_adj * & + (1.0 + sin(solar / 2.0)**2) * abs(cos(aspect - sazi)) +! sin(0.5*theta)**2 = 0.5 + 0.5*cos(theta) !--------------------------------------------------------------- ! sky diffuse irradiation for anisotropical surface ! see Iqbal, 1983 "an introduction to solar ! radiation" Hay model - edif_t = edif_h(i, j) * (ts(i, j) * cos(it) / cos(solar) + & - vd * (1 - ts(i, j))) + eadj - rdir = edir_t / (edir_h(i, j) + edif_h(i, j)) - rdif = edif_t / (edir_h(i, j) + edif_h(i, j)) - rtotal = (edir_t + edif_t) / (edir_h(i, j) + edif_h(i, j)) - - rth = (rori - s_mod(i, j) * ref_lm(i)) / (1 - s_mod(i, j) * & - ref_lm(i)) - - if (rtotal .le. rth) then - bb_angle = fs(i, j) / cos(solar) + (1 - fs(i, j)) * & - ts(i, j) / cos(solar) - cc_angle = -rth + (1.0 - fs(i, j)) * vd * & - (1.0 - ts(i, j)) + & - eadj / (edir_h(i, j) + edif_h(i, j)) - tttt = -cc_angle / bb_angle - if (tttt .gt. 1.0) tttt = 1.0 - if (tttt .lt. -1.0) tttt = -1.0 - angle_th = acos(tttt) * 180.0 / pi - angle = 90.0 - it_angle(i, j) + angle_th - edir_t = edir_h(i, j) * (cos(it) + cos(angle * pib)) / & - (cos(solar) +cos(angle * pib)) - rdir = edir_t/(edir_h(i, j)+edif_h(i, j)) - rtotal = (edir_t+edif_t)/(edir_h(i, j)+edif_h(i, j)) - endif + edif_t = tile_edif_h * (tile_ts * cosit / cossolar + & + vd * (1 - tile_ts)) + eadj + rdir = edir_t / (tile_edir_h + tile_edif_h) + rdif = edif_t / (tile_edir_h + tile_edif_h) + rtotal = (edir_t + edif_t) / (tile_edir_h + tile_edif_h) + + rth = (rori - tile_s_mod * tile_ref_lm) / (1 - tile_s_mod * & + tile_ref_lm) + + +! if (rtotal .le. rth) then + bb_angle = tile_fs / cossolar + (1 - tile_fs) * & + tile_ts / cossolar + cc_angle = -rth + (1.0 - tile_fs) * vd * & + (1.0 - tile_ts) + & + eadj / (tile_edir_h + tile_edif_h) +! tttt = -cc_angle / bb_angle + tttt = max(min(-cc_angle / bb_angle,1.0),-1.0) +! if (tttt .gt. 1.0) tttt = 1.0 +! if (tttt .lt. -1.0) tttt = -1.0 +! angle_th = acos(tttt) * 180.0 / pi + angle_th = acos(tttt) +!!! Do this in radians instead of degrees +! angle = 90.0 - it_angle(i, j) + angle_th + angle = pi / 2.0 - it + angle_th +! edir_t = edir_h(i, j) * (cos(it) + cos(angle * pib)) / & +! (cos(solar) +cos(angle * pib)) + edir_t = tile_edir_h * ( cosit + cos(angle) ) / & + ( cossolar + cos(angle) ) + rdir_tmp = edir_t/(tile_edir_h+tile_edif_h) + rtotal_tmp = (edir_t+edif_t)/(tile_edir_h+tile_edif_h) + where( rtotal .le. rth ) + rtotal = rtotal_tmp + rdir = rdir_tmp + end where +! endif !---------------------------------------------------------------- ! calculate normalized BRDF shape function for sloping surface - ann_s = RL_brdf(it_brdf,et_brdf,ra_sl,hb,br,1.0,norm_1,norm_2) + ann_s = RL_brdf(it_brdf,et_brdf,ra_sl,hb,br,1.0,norm_1,norm_2,veclen) + !---------------------------------------------------------------- ! calculate black sky albedo for sloar angle - aa_solars = black_sky(1.0,norm_1,norm_2,it_bk) + aa_solars = black_sky(1.0,norm_1,norm_2,it_bk,veclen) !-------------------------------------------------------------- ! calculate black sky albedo for view angle - aa_views = black_sky(1.0,norm_1,norm_2,et_bk) + aa_views = black_sky(1.0,norm_1,norm_2,et_bk,veclen) !------------------------------------------------------------- aa_slope = & - (rdir * (fv(i, j) * ann_s + (1.0-fv(i, j)) * & - aa_solars) + rdif * (fv(i, j) * aa_views + & - (1.0-fv(i, j)) * aa_white)) / aa_white - - a_eqs = (rtotal - aa_slope) * s_mod(i, j) * & - (1 - s_mod(i, j) * ref_lm(i)) - b_eqs = aa_slope + ref_lm(i) * (1 - aa_slope) * s_mod(i, j) - c_eqs = -ref_lm(i) - - if (abs(a_eqs) .lt. 0.00001) then - ref_bars = -c_eqs / b_eqs - else - ref_bars = (-b_eqs + sqrt(b_eqs**2 -4 * a_eqs * c_eqs)) / & - (2*a_eqs) - endif + (rdir * (tile_fv * ann_s + (1.0-tile_fv) * & + aa_solars) + rdif * (tile_fv * aa_views + & + (1.0-tile_fv) * aa_white)) / aa_white + + a_eqs = (rtotal - aa_slope) * tile_s_mod * & + (1 - tile_s_mod * tile_ref_lm) + b_eqs = aa_slope + tile_ref_lm * (1 - aa_slope) * tile_s_mod +! c_eqs = -ref_lm(i) + +! if (abs(a_eqs) .lt. 0.00001) then +! ref_bars = -c_eqs / b_eqs +! else + ref_bars = (-b_eqs + sqrt(b_eqs**2 +4 * a_eqs * tile_ref_lm)) / & + (2*a_eqs) + where( abs(a_eqs) .lt. 0.00001 ) ref_bars = tile_ref_lm / b_eqs +! endif ref_brdfreals = ann_s * ref_bars / aa_white - ref_terrain(i) = ref_bars * fnn / aa_white - iref_terrain(i, j) = int(ref_terrain(i) * 10000.0 + 0.5) - - if (ref_terrain(i) .ge. 1) then - ref_terrain(i) = 1.0 - iref_terrain(i, j) = int(ref_terrain(i) * 10000.0 + 0.5) - endif + +! ref_terrain(i) = ref_bars * fnn / aa_white +! iref_terrain(i, j) = int(ref_terrain(i) * 10000.0 + 0.5) + tile_ref_terrain = min(ref_bars * fnn / aa_white, 1.0) + tile_iref_terrain = nint(tile_ref_terrain * 10000.0,2) +! if (ref_terrain(i) .ge. 1) then +! ref_terrain(i) = 1.0 +! iref_terrain(i, j) = int(ref_terrain(i) * 10000.0 + 0.5) +! endif ! Should test for these cases in initial tests! (ie must be lt these) ! presently comments as test for ge 90 in initial one ! - if (it_angle(i, j) .ge. 85.0) then - iref_terrain(i, j) = -999 - endif - - if (et_angle(i, j) .ge. 85.0) then - iref_terrain(i, j) = -999 - endif +! if (it_angle(i, j) .ge. 85.0) then +! iref_terrain(i, j) = -999 +! endif +! +! if (et_angle(i, j) .ge. 85.0) then +! iref_terrain(i, j) = -999 +! endif ! if in masked or otherwise unused area you get here ! put in -999 values to indicate null data - else - iref_terrain(i, j) = -999 - endif - else - ref_lm(i) = -999.0 - ref_brdf(i) = -999.0 - ref_terrain(i) = -999.0 - iref_lm(i, j) = -999 - iref_brdf(i, j) = -999 - iref_terrain(i, j) = -999 - endif - enddo +! else +! iref_terrain(i, j) = -999 +! endif + where( it .ge. 85.0*pib .OR. & + et .ge. 85.0*pib .OR. & + tile_mask_self .le. 0 .OR. & + tile_mask_castview .le. 0 .OR. & + tile_mask_castsun .le. 0 ) & + tile_iref_terrain = -999 + else + tile_iref_terrain = -999 + end if + + where (tile_a_mod(1:tilelen) .ge. 0 .and. tile_dn_1(1:tilelen) .gt. 0 ) + iref_lm(istart:iend,j) = tile_iref_lm(1:tilelen) + iref_brdf(istart:iend,j) = tile_iref_brdf(1:tilelen) + iref_terrain(istart:iend,j) = tile_iref_terrain(1:tilelen) + elsewhere +! ref_lm(i) = -999.0 +! ref_brdf(i) = -999.0 +! ref_terrain(i) = -999.0 + iref_lm(istart:iend, j) = -999 + iref_brdf(istart:iend, j) = -999 + iref_terrain(istart:iend, j) = -999 + end where + else + iref_lm(istart:iend, j) = -999 + iref_brdf(istart:iend, j) = -999 + iref_terrain(istart:iend, j) = -999 + end if + enddo enddo -END SUBROUTINE reflectance + END SUBROUTINE reflectance diff --git a/src/Makefile b/src/Makefile index 80efaeb..e94ba8b 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,7 +1,7 @@ SHELL := /bin/bash -CC := gcc -FC := gfortran -FCFLAGS := -Wall -O3 -fbounds-check +CC := icc +FC := ifort +FCFLAGS := -O3 -check bounds -assume byterecl -xHost -ipo CFLAGS := -O3 EXES := calculate_coefficients \ diff --git a/workflow/nbar.cfg b/workflow/nbar.cfg index ca8ebf4..f480662 100644 --- a/workflow/nbar.cfg +++ b/workflow/nbar.cfg @@ -48,7 +48,7 @@ tc_intermediates = TC_Intermediates rfl_output_dir = Reflectance_Outputs header_slope_target = SLOPE_ANGLE_INPUTS x_tile_size = -1 -y_tile_size = 100 +y_tile_size = 1120 [brdf] factor = 1.0 From 45cf4133269e512fb0bdb8fb68ac3680f18be273 Mon Sep 17 00:00:00 2001 From: Dale Roberts Date: Mon, 6 Jun 2016 11:46:55 +1000 Subject: [PATCH 3/5] Vector version of surface_reflectance.f90 (and called functions) discussed in meeting on Friday 3/6/16. Also included are changes in the makefiles necessary for an Intel compiler build, and updated nbar.cfg with optimal y_tile size. --- gaip/reflectance_mod.f90 | 99 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 99 insertions(+) create mode 100644 gaip/reflectance_mod.f90 diff --git a/gaip/reflectance_mod.f90 b/gaip/reflectance_mod.f90 new file mode 100644 index 0000000..d04bf5a --- /dev/null +++ b/gaip/reflectance_mod.f90 @@ -0,0 +1,99 @@ +module reflectance_mod + + real, parameter :: pi = 3.14159265358979 + real, parameter :: pib = 0.017453292519943 ! pi / 180 + +contains + + function RL_brdf (solar,view,ra,hb,br,brdf0,brdf1,brdf2, veclen) +! calculate the normalised BRDF shape function + integer :: veclen + real solar(veclen), view(veclen), ra(veclen), hb, br + real brdf0,brdf1,brdf2 + real RL_brdf(veclen) + + real rs_thick(veclen),li_sparse(veclen),secsolar(veclen),secvia(veclen) + real cossolar(veclen),cosvia(veclen),cosra(veclen),sinsolar(veclen),sinvia(veclen),sinra(veclen) + real cosxi(veclen),xi(veclen),tansolar(veclen),tanvia(veclen),theta_new_v(veclen),theta_new_s(veclen) + real d_li2(veclen),x_li(veclen),cosl(veclen),l_li(veclen),o_li(veclen) + + real tan_theta_new_v(veclen), cos_theta_new_v(veclen), sin_theta_new_v(veclen) + real tan_theta_new_s(veclen), cos_theta_new_s(veclen), sin_theta_new_s(veclen) + + RL_brdf=0.0 + +! pi = 4 * atan(1.0) +! calculate Ross-thick kernel + cossolar = cos(solar) + cosvia = cos(view) + cosra = cos(ra) + sinsolar = sin(solar) + sinvia = sin(view) + sinra = sin(ra) + cosxi = min(cossolar*cosvia+sinsolar*sinvia*cosra,1.0) +! if (cosxi .ge.1) then +! cosxi = 1 +! endif + xi = acos(cosxi) + rs_thick = ((pi/2.0-xi)*cosxi+sin(xi)) / (cossolar+cosvia) - pi/4.0 + +! alculate Li-sparse + tansolar = sinsolar/cossolar + tanvia = sinvia/cosvia +! theta_new_v = atan(br*tanvia) +! theta_new_s = atan(br*tansolar) + tan_theta_new_v = br*tanvia + ! Trig identities + cos_theta_new_v = 1.0 / sqrt(1.0 + tan_theta_new_v * tan_theta_new_v) + sin_theta_new_v = tan_theta_new_v * cos_theta_new_v + + tan_theta_new_s = br*tansolar + ! Trig identities + cos_theta_new_s = 1.0 / sqrt(1.0 + tan_theta_new_s * tan_theta_new_s) + sin_theta_new_s = tan_theta_new_s * cos_theta_new_s + +! cosxi = cos(theta_new_s)*cos(theta_new_v)+sin(theta_new_s)*sin(theta_new_v)*cosra + cosxi = min( cos_theta_new_s*cos_theta_new_v + sin_theta_new_s*sin_theta_new_v*cosra, 1.0 ) +! if (cosxi.ge.1) then +! cosxi = 1 +! endif + secsolar = 1.0/cos_theta_new_s + secvia = 1.0/cos_theta_new_v +! d_li2 = abs(tan(theta_new_s)**2+tan(theta_new_v)**2 - 2.0*tan(theta_new_s)*tan(theta_new_v)*cosra) + d_li2 = abs(tan_theta_new_s**2+tan_theta_new_v**2 - 2.0*tan_theta_new_s*tan_theta_new_v*cosra) + x_li = tan_theta_new_s*tan_theta_new_v*sinra + cosl = hb*sqrt(d_li2+x_li**2)/(secsolar+secvia) +! if (cosl .ge. 1.0) then +! o_li=0.0 +! else + l_li=acos(cosl) + o_li=(l_li-sin(l_li)*cosl)*(secsolar+secvia)/pi + where( cosl .ge. 1.0 ) o_li = 0.0 +! endif + li_sparse=o_li-(secsolar+secvia)+0.5*(1.0+cosxi) * secsolar*secvia + RL_brdf=brdf0+brdf1*rs_thick+brdf2*li_sparse + return + end function RL_brdf + + function black_sky(brdf0,brdf1,brdf2,theta,veclen) + integer veclen + real brdf0,brdf1,brdf2 + real theta(veclen) + real black_sky(veclen) + +! theta should be in radians + black_sky = brdf0 + & + brdf1*(-0.007574-0.070987*theta**2 + 0.307588*theta**3) + & + brdf2*(-1.284909-0.166314*theta**2 + 0.041840*theta**3) + return + end function black_sky + + + function white_sky(brdf0, brdf1, brdf2) + real brdf0, brdf1, brdf2 + real white_sky + white_sky = brdf0 + brdf1*0.189184 - brdf2*1.377622 + return + end function white_sky + +end module reflectance_mod From d02add9d79c341de67a1a407926ccbe153c9c8d5 Mon Sep 17 00:00:00 2001 From: Dale Roberts Date: Mon, 6 Jun 2016 12:13:05 +1000 Subject: [PATCH 4/5] Updated working copy to latest revision, merged changes --- .gitignore | 1 + gaip/Makefile | 4 +- gaip/black_sky.f90 | 6 +- gaip/brdf_shape.f90 | 82 ++++-- gaip/setup.py | 4 +- gaip/surface_reflectance.f90 | 537 +++++++++++++++++++++++------------ src/Makefile | 6 +- workflow/nbar.cfg | 2 +- 8 files changed, 414 insertions(+), 228 deletions(-) diff --git a/.gitignore b/.gitignore index 6df902f..b72b97c 100644 --- a/.gitignore +++ b/.gitignore @@ -35,3 +35,4 @@ nosetests.xml .idea *.iml *.ipr +*~ \ No newline at end of file diff --git a/gaip/Makefile b/gaip/Makefile index 3be082d..be5723d 100644 --- a/gaip/Makefile +++ b/gaip/Makefile @@ -1,4 +1,4 @@ -FC := gfortran +FC := ifort F2PY:= f2py all: _cast_shadow_mask.so _slope_self_shadow.so _surface_reflectance.so sys_variables.o sys_variables.mod _sat_sol_angles.so _satellite_model.so _track_time_info.so unmiximage.pyf unmiximage.so _interpolation.so @@ -6,7 +6,7 @@ all: _cast_shadow_mask.so _slope_self_shadow.so _surface_reflectance.so sys_vari _cast_shadow_mask.so _slope_self_shadow.so _surface_reflectance.so sys_variables.o sys_variables.mod _sat_sol_angles.so _satellite_model.so _track_time_info.so: $$(module load gcc/5.2.0) $(FC) -c sys_variables.f90 - python setup.py build_ext -i + python setup.py config_fc --opt "-O3 -g" --arch "-xHost" build_ext -i unmiximage.pyf unmiximage.so: $(F2PY) -h unmiximage.pyf -m unmiximage unmiximage.f90 diff --git a/gaip/black_sky.f90 b/gaip/black_sky.f90 index 4168623..40f3bbe 100644 --- a/gaip/black_sky.f90 +++ b/gaip/black_sky.f90 @@ -1,6 +1,8 @@ -real function black_sky(brdf0,brdf1,brdf2,theta) +function black_sky(brdf0,brdf1,brdf2,theta,veclen) + integer veclen real brdf0,brdf1,brdf2 - real theta + real theta(veclen) + real black_sky(veclen) ! theta should be in radians black_sky = brdf0 + & diff --git a/gaip/brdf_shape.f90 b/gaip/brdf_shape.f90 index 76aeba0..a36c33a 100644 --- a/gaip/brdf_shape.f90 +++ b/gaip/brdf_shape.f90 @@ -1,16 +1,23 @@ -real function RL_brdf (solar,view,ra,hb,br,brdf0,brdf1,brdf2) +function RL_brdf (solar,view,ra,hb,br,brdf0,brdf1,brdf2,veclen) ! calculate the normalised BRDF shape function - real pi - real solar, view, ra, hb, br + integer veclen + real solar(veclen), view(veclen), ra(veclen), hb, br real brdf0,brdf1,brdf2 - real rs_thick,li_sparse,secsolar,secvia - real cossolar,cosvia,cosra,sinsolar,sinvia,sinra - real cosxi,xi,tansolar,tanvia,theta_new_v,theta_new_s - real d_li2,x_li,cosl,l_li,o_li + real RL_brdf(veclen) + + real rs_thick(veclen),li_sparse(veclen),secsolar(veclen),secvia(veclen) + real cossolar(veclen),cosvia(veclen),cosra(veclen),sinsolar(veclen),sinvia(veclen),sinra(veclen) + real cosxi(veclen),xi(veclen),tansolar(veclen),tanvia(veclen),theta_new_v(veclen),theta_new_s(veclen) + real d_li2(veclen),x_li(veclen),cosl(veclen),l_li(veclen),o_li(veclen) + + real tan_theta_new_v(veclen), cos_theta_new_v(veclen), sin_theta_new_v(veclen) + real tan_theta_new_s(veclen), cos_theta_new_s(veclen), sin_theta_new_s(veclen) + + real, parameter :: pi = 3.14159265358979 RL_brdf=0.0 - pi = 4 * atan(1.0) +! pi = 4 * atan(1.0) ! calculate Ross-thick kernel cossolar = cos(solar) cosvia = cos(view) @@ -18,34 +25,47 @@ real function RL_brdf (solar,view,ra,hb,br,brdf0,brdf1,brdf2) sinsolar = sin(solar) sinvia = sin(view) sinra = sin(ra) - cosxi = cossolar*cosvia+sinsolar*sinvia*cosra - if (cosxi .ge.1) then - cosxi = 1 - endif + cosxi = min(cossolar*cosvia+sinsolar*sinvia*cosra,1.0) +! if (cosxi .ge.1) then +! cosxi = 1 +! endif xi = acos(cosxi) - rs_thick = ((pi/2.0-xi)*cos(xi)+sin(xi)) / (cossolar+cosvia) - pi/4.0 + rs_thick = ((pi/2.0-xi)*cosxi+sin(xi)) / (cossolar+cosvia) - pi/4.0 ! alculate Li-sparse tansolar = sinsolar/cossolar tanvia = sinvia/cosvia - theta_new_v = atan(br*tanvia) - theta_new_s = atan(br*tansolar) - cosxi = cos(theta_new_s)*cos(theta_new_v)+sin(theta_new_s)*sin(theta_new_v)*cosra - if (cosxi.ge.1) then - cosxi = 1 - endif - secsolar = 1.0/cos(theta_new_s) - secvia = 1.0/cos(theta_new_v) - d_li2 = abs(tan(theta_new_s)**2+tan(theta_new_v)**2 - 2.0*tan(theta_new_s)*tan(theta_new_v)*cosra) - x_li = tan(theta_new_s)*tan(theta_new_v)*sinra +! theta_new_v = atan(br*tanvia) +! theta_new_s = atan(br*tansolar) + tan_theta_new_v = br*tanvia + ! Trig identities + cos_theta_new_v = 1.0 / sqrt(1 + tan_theta_new_v * tan_theta_new_v) + sin_theta_new_v = tan_theta_new_v * cos_theta_new_v + + tan_theta_new_s = br*tansolar + ! Trig identities + cos_theta_new_s = 1.0 / sqrt(1 + tan_theta_new_s * tan_theta_new_s) + sin_theta_new_s = tan_theta_new_s * cos_theta_new_s + +! cosxi = cos(theta_new_s)*cos(theta_new_v)+sin(theta_new_s)*sin(theta_new_v)*cosra + cosxi = min( cos_theta_new_s*cos_theta_new_v + sin_theta_new_s*sin_theta_new_s*cosra, 1.0 ) +! if (cosxi.ge.1) then +! cosxi = 1 +! endif + secsolar = 1.0/cos_theta_new_s + secvia = 1.0/cos_theta_new_v +! d_li2 = abs(tan(theta_new_s)**2+tan(theta_new_v)**2 - 2.0*tan(theta_new_s)*tan(theta_new_v)*cosra) + d_li2 = abs(tan_theta_new_s**2+tan_theta_new_v**2 - 2.0*tan_theta_new_s*tan_theta_new_v*cosra) + x_li = tan_theta_new_s*tan_theta_new_v*sinra cosl = hb*sqrt(d_li2+x_li**2)/(secsolar+secvia) - if (cosl .ge. 1.0) then - o_li=0.0 - else - l_li=acos(cosl) - o_li=(l_li-sin(l_li)*cos(l_li))*(secsolar+secvia)/pi - endif - li_sparse=o_li-(secsolar+secvia)+0.5*(1.0+cosxi) * secsolar*secvia - RL_brdf=brdf0+brdf1*rs_thick+brdf2*li_sparse +! if (cosl .ge. 1.0) then +! o_li=0.0 +! else + l_li=acos(cosl) + o_li=(l_li-sin(l_li)*cosl)*(secsolar+secvia)/pi + where( cosl .ge. 1.0 ) o_li = 0.0 +! endif + li_sparse=o_li-(secsolar+secvia)+0.5*(1.0+cosxi) * secsolar*secvia + RL_brdf=brdf0+brdf1*rs_thick+brdf2*li_sparse return end diff --git a/gaip/setup.py b/gaip/setup.py index 4501329..782aaf7 100644 --- a/gaip/setup.py +++ b/gaip/setup.py @@ -25,9 +25,7 @@ 'geo2metres_pixel_size.f90']), Extension(name='_surface_reflectance', sources=['surface_reflectance.f90', - 'white_sky.f90', - 'black_sky.f90', - 'brdf_shape.f90']), + 'reflectance_mod.f90']), Extension(name='_satellite_model', sources=['geo2metres_pixel_size.f90', 'satellite_model.f90']), diff --git a/gaip/surface_reflectance.f90 b/gaip/surface_reflectance.f90 index 4b14793..2b6e0bc 100644 --- a/gaip/surface_reflectance.f90 +++ b/gaip/surface_reflectance.f90 @@ -36,7 +36,9 @@ SUBROUTINE reflectance( & ! Calculates lambertian, brdf corrected and terrain corrected surface ! reflectance. + use reflectance_mod + implicit none ! input parameters integer nrow, ncol ! we should be passing in transposed arrays cols = rows and vice versa real*4 rori ! threshold for terrain correction @@ -89,29 +91,61 @@ SUBROUTINE reflectance( & !f2py intent(inout) iref_lm, iref_brdf, iref_terrain ! internal parameters + real, parameter :: hb = 2.0 + real, parameter :: br = 1.0 + integer, parameter :: veclen = 16 + + !!! Tiled input quantities + real*4 tile_dn_1(veclen) ! raw image (in real4 as replacement for dn array) + integer*2 tile_mask_self(veclen) ! mask + integer*2 tile_mask_castsun(veclen) ! self shadow mask + integer*2 tile_mask_castview(veclen) ! cast shadow mask +! real*4 tile_solar_angle(veclen) ! solar zenith angle +! real*4 tile_sazi_angle(veclen) ! solar azimuth angle +! real*4 tile_view_angle(veclen) ! view angle (for flat surface) + real*4 tile_rela_angle(veclen) ! relative azimuth angle (for flat surface) +! real*4 tile_slope_angle(veclen) ! slop angle +! real*4 tile_aspect_angle(veclen) ! aspect angle +! real*4 tile_it_angle(veclen) ! incident angle (for inclined surface) +! real*4 tile_et_angle(veclen) ! exiting angle (for inclined surface) + real*4 tile_rela_slope(veclen) ! relative angle (for inclined surface) + real*4 tile_a_mod(veclen) ! modtran output (a) + real*4 tile_b_mod(veclen) ! modtran output (b) + real*4 tile_s_mod(veclen) ! modtran output (s) + real*4 tile_fv(veclen) ! modtran output (fs) + real*4 tile_fs(veclen) ! modtran output (fv) + real*4 tile_ts(veclen) ! modtran output (ts) + real*4 tile_edir_h(veclen) ! modtran output (direct irradiance) + real*4 tile_edif_h(veclen) ! modtran output (diffuse irradiance) + integer*2 tile_iref_lm(veclen) ! atmospheric corrected lambertial reflectance + integer*2 tile_iref_brdf(veclen) ! atmospheric and brdf corrected reflectance + integer*2 tile_iref_terrain(veclen) ! atmospheric and brdf and terrain correc + + real*4 tile_ref_lm(veclen) + real*4 tile_ref_brdf(veclen) + real*4 tile_ref_terrain(veclen) +!!! Tile sizing + integer :: istart, iend, tilelen + integer i, j - real pi, pib - real ann_f, aa_viewf, aa_solarf, aa_white - real ann_s, aa_views, aa_solars - real lt, norm_1, norm_2 - real solar, sazi, view, slope, aspect, ra_lm, ra_sl, it, et - - double precision a_eqf, b_eqf, c_eqf, ref_barf, aa_flat - double precision a_eqs, b_eqs, c_eqs, ref_bars, aa_slope - double precision ref_brdfrealf, ref_brdfreals - double precision bb_angle, cc_angle, tttt - real hb, br, fnn, it_brdf, it_bk, et_brdf, et_bk - real angle_th - real vt, vd,angle, edir_t, eadj, edif_t, rdir, rdif, rtotal - real RL_brdf, black_sky, white_sky - real cosslope, rth - external RL_brdf, black_sky, white_sky + real ann_f(veclen), aa_viewf(veclen), aa_solarf(veclen), aa_white + real ann_s(veclen), aa_views(veclen), aa_solars(veclen) + real lt(veclen), norm_1, norm_2 + real solar(veclen), sazi(veclen), view(veclen), slope(veclen), aspect(veclen), ra_lm(veclen), ra_sl(veclen), it(veclen), et(veclen) + + double precision a_eqf(veclen), b_eqf(veclen), c_eqf, ref_barf(veclen), aa_flat(veclen) + double precision a_eqs(veclen), b_eqs(veclen), c_eqs, ref_bars(veclen), aa_slope(veclen) + double precision ref_brdfrealf(veclen), ref_brdfreals(veclen) + double precision bb_angle(veclen), cc_angle(veclen), tttt(veclen) + real fnn, fnn_a(1), it_brdf(veclen), it_bk(veclen), et_brdf(veclen), et_bk(veclen) + real angle_th(veclen) + real vt(veclen), vd(veclen),angle(veclen), edir_t(veclen), eadj(veclen), edif_t(veclen), rdir(veclen), rdif(veclen), rtotal(veclen), rtotal_tmp(veclen), rdir_tmp(veclen) + real cossolar(veclen),cosit(veclen),rth(veclen) + + logical :: tilemask(veclen) ! li-sparse parameters - hb = 2.0 - br = 1.0 - pi = 4 * atan(1.0) - pib = pi / 180.0 + norm_1 = brdf1 / brdf0 norm_2 = brdf2 / brdf0 @@ -119,228 +153,359 @@ SUBROUTINE reflectance( & ! calculate white sky albedo aa_white = white_sky(1.0, norm_1, norm_2) ! calcualte BRDF at 45 solar angle and 0 view angle - fnn = RL_brdf(45 * pib, 0.0, 0.0, hb, br, 1.0, norm_1, norm_2) -! print*,fnn + fnn_a = RL_brdf( (/ 45 * pib /), (/ 0.0 /), (/ 0.0 /), hb, br, 1.0, norm_1, norm_2,1) + fnn = fnn_a(1) ! Now loop over the cols of the images do j=1,ncol !---------------------------------------------------------------------- ! convert byte to integer for the raw image - do i=1,nrow - if (dn_1(i, j) .lt. 0) then - dn(i) = dn_1(i, j) + 65536 - else - dn(i) = dn_1(i, j) - endif - enddo +! do i=1,nrow +! if (dn_1(i, j) .lt. 0) then +! dn(i) = dn_1(i, j) + 65536 +! else +! dn(i) = dn_1(i, j) +! endif +! enddo + istart = 0 + iend = 0 + do while ( iend < nrow ) + istart = iend + 1 + iend = iend + veclen + if ( iend > nrow ) iend = nrow + tilelen = iend - istart + 1 + + where ( dn_1(istart:iend,j) .lt. 0 ) + tile_dn_1(1:tilelen) = real(dn_1(istart:iend,j) + 65536,4) + elsewhere + tile_dn_1(1:tilelen) = real(dn_1(istart:iend,j),4) + end where + + + if( ANY( a_mod(istart:iend,j) .gt. 0 .AND. tile_dn_1(1:tilelen) .gt. 0.0 ) ) then + + tile_a_mod(1:tilelen) = a_mod(istart:iend,j) + tile_b_mod(1:tilelen) = b_mod(istart:iend,j) + tile_s_mod(1:tilelen) = s_mod(istart:iend,j) + + tile_fs(1:tilelen) = fs(istart:iend,j) + tile_fv(1:tilelen) = fv(istart:iend,j) + tile_ts(1:tilelen) = ts(istart:iend,j) + + tile_edir_h(1:tilelen) = edir_h(istart:iend,j) + tile_edif_h(1:tilelen) = edif_h(istart:iend,j) +!!! Mask tiles + tile_mask_self(1:tilelen) = mask_self(istart:iend,j) + tile_mask_castsun(1:tilelen) = mask_castsun(istart:iend,j) + tile_mask_castview(1:tilelen) = mask_castview(istart:iend,j) + ! now loop over the rows of the images - do i=1,nrow +! do i=1,nrow ! if valid masks and valid digital number then do the calcs - if (a_mod(i, j) .ge. 0 .and. dn(i) .gt. 0) then - if (rela_angle(i, j) .gt. 180) rela_angle(i, j) = rela_angle(i, j) - 360 - if (rela_slope(i, j) .gt. 180) rela_slope(i, j) = rela_slope(i, j) - 360 - +!!! Mask after calculations +! if (a_mod(i, j) .ge. 0 .and. dn(i) .gt. 0) then + +! if (rela_angle(i, j) .gt. 180) rela_angle(i, j) = rela_angle(i, j) - 360 +! if (rela_slope(i, j) .gt. 180) rela_slope(i, j) = rela_slope(i, j) - 360 + where ( rela_angle(istart:iend,j) .gt. 180 ) + tile_rela_angle(1:tilelen) = rela_angle(istart:iend,j) - 360 + elsewhere + tile_rela_angle(1:tilelen) = rela_angle(istart:iend,j) + end where + + where ( rela_slope(istart:iend,j) .gt. 180 ) + tile_rela_slope(1:tilelen) = rela_slope(istart:iend,j) - 360 + elsewhere + tile_rela_slope(1:tilelen) = rela_slope(istart:iend,j) + end where + ! convert angle to radians - solar = solar_angle(i, j) * pib - sazi = sazi_angle(i, j) * pib - view = view_angle(i, j) * pib - slope = slope_angle(i, j) * pib - aspect = aspect_angle(i, j) * pib - ra_lm = rela_angle(i, j) * pib - ra_sl = rela_slope(i, j) * pib - it = it_angle(i, j) * pib - et = et_angle(i, j) * pib - - if (it_angle(i, j) .ge. 70.0) then - it_brdf= 70.0 * pib - else - it_brdf = it_angle(i, j) * pib - endif - - if (it_angle(i, j) .ge. 80.0) then - it_bk = 80.0 * pib - else - it_bk = it_angle(i, j) * pib - endif - - if (et_angle(i, j) .ge. 60.0) then - et_brdf = 60.0 * pib - else - et_brdf = et_angle(i, j) * pib - endif - - if (et_angle(i, j) .ge. 80.0) then - et_bk = 80.0 * pib - else - et_bk = et_angle(i, j) * pib - endif +! solar = solar_angle(i, j) * pib + solar(1:tilelen) = solar_angle(istart:iend,j) + solar = solar * pib + ! Trig is expensive + cossolar = cos(solar) +! sazi = sazi_angle(i, j) * pib + sazi(1:tilelen) = sazi_angle(istart:iend,j) + sazi = sazi * pib + ! view = view_angle(i, j) * pib + view(1:tilelen) = view_angle(istart:iend,j) + view = view * pib + ! slope = slope_angle(i, j) * pib + slope(1:tilelen) = slope_angle(istart:iend,j) + slope = slope * pib + ! aspect = aspect_angle(i, j) * pib + aspect(1:tilelen) = aspect_angle(istart:iend,j) + aspect = aspect * pib + ! ra_lm = rela_angle(i, j) * pib + ra_lm(1:tilelen) = rela_angle(istart:iend,j) + ra_lm = ra_lm * pib + ! ra_sl = rela_slope(i, j) * pib + ra_sl(1:tilelen) = rela_slope(istart:iend,j) + ra_sl = ra_sl * pib + ! it = it_angle(i, j) * pib + it(1:tilelen) = it_angle(istart:iend,j) + ! et = et_angle(i, j) * pib + et(1:tilelen) = et_angle(istart:iend,j) + +! if (it_angle(i, j) .ge. 70.0) then +! it_brdf= 70.0 * pib +! else +! it_brdf = it_angle(i, j) * pib +! endif + it_brdf = min(70.0,it) * pib + +! if (it_angle(i, j) .ge. 80.0) then +! it_bk = 80.0 * pib +! else +! it_bk = it_angle(i, j) * pib +! endif + it_bk = min(80.0, it) * pib + + it = it * pib + ! Trig is expensive + cosit = cos(it) + + ! if (et_angle(i, j) .ge. 60.0) then + ! et_brdf = 60.0 * pib + ! else + ! et_brdf = et_angle(i, j) * pib + ! endif + et_brdf = min(60.0,et) * pib + +! if (et_angle(i, j) .ge. 80.0) then +! et_bk = 80.0 * pib +! else +! et_bk = et_angle(i, j) * pib +! endif + et_bk = min(80.0,et) * pib + + et = et * pib !------------------------------------------------ ! for flat surface ! calculate radiance at top atmosphere - lt = bias + dn(i) * slope_ca +! lt = bias + dn(i) * slope_ca + lt = bias + tile_dn_1 * slope_ca ! calcualte lambetian reflectance with bilinear average - ref_lm(i) = (lt - b_mod(i, j)) / (a_mod(i, j) + s_mod(i, j) * & - (lt - b_mod(i, j))) - iref_lm(i, j) = ref_lm(i) * 10000 + 0.5 +! ref_lm(i) = (lt - b_mod(i, j)) / (a_mod(i, j) + s_mod(i, j) * & +! (lt - b_mod(i, j))) +! iref_lm(i, j) = ref_lm(i) * 10000 + 0.5 + tile_ref_lm = (lt - tile_b_mod) / ( tile_a_mod + tile_s_mod * & + ( lt - tile_b_mod )) + tile_iref_lm = nint(tile_ref_lm * 10000.0,2) ! set as small number if atmospheric corrected reflectance ! below 0.001 - if (ref_lm(i).lt. 0.001) then - ref_lm(i) = 0.001 - ref_brdf(i) = 0.001 - iref_lm(i, j) = 10 - iref_brdf(i, j) = 10 - else +!!! Set these at the end +! if (ref_lm(i).lt. 0.001) then +! ref_lm(i) = 0.001 +! ref_brdf(i) = 0.001 +! iref_lm(i, j) = 10 +! iref_brdf(i, j) = 10 + ! else ! calculate normalized BRDF shape function - ann_f = RL_brdf(solar, view, ra_lm, hb, br, 1.0, norm_1, norm_2) + if( all(tile_ref_lm(1:tilelen) .lt.0.001) ) then + tile_ref_lm = 0.001 + tile_ref_brdf = 0.001 + tile_iref_lm = 10 + tile_iref_brdf = 10 + else + ann_f = RL_brdf(solar, view, ra_lm, hb, br, 1.0, norm_1, norm_2,veclen) ! calculate black sky albedo for sloar angle - aa_solarf = black_sky(1.0, norm_1, norm_2, solar) + aa_solarf = black_sky(1.0, norm_1, norm_2, solar,veclen) ! calculate black sky albedo for view angle - aa_viewf = black_sky(1.0, norm_1, norm_2, view) + aa_viewf = black_sky(1.0, norm_1, norm_2, view,veclen) ! - aa_flat = (fv(i, j) * (fs(i, j) * ann_f + (1.0 - fs(i, j)) * & - aa_viewf) + (1.0 - fv(i, j)) * (fs(i, j) * & - aa_solarf + (1.0 - fs(i, j)) * aa_white)) / aa_white - - a_eqf = (1 - aa_flat) * s_mod(i, j) * (1 - s_mod(i, j) * & - ref_lm(i)) - b_eqf = aa_flat + ref_lm(i) * (1 - aa_flat) * s_mod(i, j) - c_eqf = -ref_lm(i) - - if (abs(a_eqf) .lt. 0.0000001) then - ref_barf = -c_eqf / b_eqf - else - ref_barf = (-b_eqf + & - sqrt(b_eqf**2 - 4 * a_eqf * c_eqf)) / & - (2 * a_eqf) - endif - ref_brdfrealf = ann_f * ref_barf / aa_white - ref_brdf(i) = ref_barf * fnn / aa_white - iref_brdf(i, j) = ref_brdf(i) * 10000 + 0.5 + aa_flat = (tile_fv * (tile_fs * ann_f + (1.0 - tile_fs) * & + aa_viewf) + (1.0 - tile_fv) * (tile_fs * & + aa_solarf + (1.0 - tile_fs) * aa_white)) / aa_white + + a_eqf = (1 - aa_flat) * tile_s_mod * (1 - tile_s_mod * & + tile_ref_lm) + b_eqf = aa_flat + tile_ref_lm * (1 - aa_flat) * tile_s_mod +! c_eqf = -ref_lm(i) + +! if (abs(a_eqf) .lt. 0.0000001) then +! ref_barf = tile_ref_lm / b_eqf +! else + ref_barf = (-b_eqf + & + sqrt(b_eqf**2 + 4 * a_eqf * tile_ref_lm)) / & + (2 * a_eqf) + where( abs(a_eqf) .lt. 0.0000001 ) ref_barf = tile_ref_lm / b_eqf +! endif + ref_brdfrealf = ann_f * ref_barf / aa_white + tile_ref_brdf = min(ref_barf * fnn / aa_white, 1.0) + tile_iref_brdf = nint(tile_ref_brdf * 10000.0,2) +! tile_iref_brdf = int(tile_ref_brdf * 10000.0 + 0.5,2) ! this is to ensure that the brdf correction -! is the same as (or as close as possible to) the original NBAR version - if (ref_brdf(i) .ge. 1) then - ref_brdf(i) = 1.0 - iref_brdf(i, j) = ref_brdf(i) * 10000 - endif - - endif + ! is the same as (or as close as possible to) the original NBAR version +!!! Condition from the top of this section. + where( tile_ref_lm .lt.0.001 ) + tile_ref_lm = 0.001 + tile_ref_brdf = 0.001 + tile_iref_lm = 10 + tile_iref_brdf = 10 + end where + end if + +! if (ref_brdf(i) .ge. 1) then +! ref_brdf(i) = 1.0 +! iref_brdf(i, j) = ref_brdf(i) * 10000 +! endif + +! endif !------------------------------------------------------------------- ! conduct terrain correction - if ((mask_self(i, j) .gt. 0) .and. & - (mask_castsun(i, j).gt. 0) .and. & - (mask_castview(i, j) .gt. 0) .and. (it_angle(i, j) .lt. 90.0) & - .and. (et_angle(i, j) .lt. 90.0)) then +!!! Tile-wide mask at the start, per-element mask at the end + tilemask(1:tilelen) = ((tile_mask_self(1:tilelen) .gt. 0) .and. & + (tile_mask_castsun(1:tilelen) .gt. 0) .and. & + (tile_mask_castview(1:tilelen) .gt. 0) .and. & + (it(1:tilelen) .lt. 85.0*pib ) .and. & + (et(1:tilelen) .lt. 85.0*pib ) ) + if ( ANY( tilemask(1:tilelen) ) ) then !---------------------------------------------------------- - cosslope = cos(slope) ! calculate vd and vt - vd = 0.5 * (1.0 + cosslope) + vd = 0.5 * (1.0 + cos(slope)) vt = 1.0 - vd !--------------------------------------------------------- ! calculate direct irradiance ! Note the account taken of threshold - - edir_t = edir_h(i, j) * cos(it) / cos(solar) + edir_t = tile_edir_h * cosit / cossolar ! calcualte adjacent irradiance for anisotropical surface ! see Iqbal, 1983 "an introduction to solar ! radiation" - eadj = (edir_h(i, j) + edif_h(i, j)) * vt * ref_adj * & - (1.0 + sin(solar / 2.0)**2) * abs(cos(aspect - sazi)) + eadj = (tile_edir_h + tile_edif_h) * vt * ref_adj * & + (1.0 + sin(solar / 2.0)**2) * abs(cos(aspect - sazi)) +! sin(0.5*theta)**2 = 0.5 + 0.5*cos(theta) !--------------------------------------------------------------- ! sky diffuse irradiation for anisotropical surface ! see Iqbal, 1983 "an introduction to solar ! radiation" Hay model - edif_t = edif_h(i, j) * (ts(i, j) * cos(it) / cos(solar) + & - vd * (1 - ts(i, j))) + eadj - rdir = edir_t / (edir_h(i, j) + edif_h(i, j)) - rdif = edif_t / (edir_h(i, j) + edif_h(i, j)) - rtotal = (edir_t + edif_t) / (edir_h(i, j) + edif_h(i, j)) - - rth = (rori - s_mod(i, j) * ref_lm(i)) / (1 - s_mod(i, j) * & - ref_lm(i)) - - if (rtotal .le. rth) then - bb_angle = fs(i, j) / cos(solar) + (1 - fs(i, j)) * & - ts(i, j) / cos(solar) - cc_angle = -rth + (1.0 - fs(i, j)) * vd * & - (1.0 - ts(i, j)) + & - eadj / (edir_h(i, j) + edif_h(i, j)) - tttt = -cc_angle / bb_angle - if (tttt .gt. 1.0) tttt = 1.0 - if (tttt .lt. -1.0) tttt = -1.0 - angle_th = acos(tttt) * 180.0 / pi - angle = 90.0 - it_angle(i, j) + angle_th - edir_t = edir_h(i, j) * (cos(it) + cos(angle * pib)) / & - (cos(solar) +cos(angle * pib)) - rdir = edir_t/(edir_h(i, j)+edif_h(i, j)) - rtotal = (edir_t+edif_t)/(edir_h(i, j)+edif_h(i, j)) - endif + edif_t = tile_edif_h * (tile_ts * cosit / cossolar + & + vd * (1 - tile_ts)) + eadj + rdir = edir_t / (tile_edir_h + tile_edif_h) + rdif = edif_t / (tile_edir_h + tile_edif_h) + rtotal = (edir_t + edif_t) / (tile_edir_h + tile_edif_h) + + rth = (rori - tile_s_mod * tile_ref_lm) / (1 - tile_s_mod * & + tile_ref_lm) + + +! if (rtotal .le. rth) then + bb_angle = tile_fs / cossolar + (1 - tile_fs) * & + tile_ts / cossolar + cc_angle = -rth + (1.0 - tile_fs) * vd * & + (1.0 - tile_ts) + & + eadj / (tile_edir_h + tile_edif_h) +! tttt = -cc_angle / bb_angle + tttt = max(min(-cc_angle / bb_angle,1.0),-1.0) +! if (tttt .gt. 1.0) tttt = 1.0 +! if (tttt .lt. -1.0) tttt = -1.0 +! angle_th = acos(tttt) * 180.0 / pi + angle_th = acos(tttt) +!!! Do this in radians instead of degrees +! angle = 90.0 - it_angle(i, j) + angle_th + angle = pi / 2.0 - it + angle_th +! edir_t = edir_h(i, j) * (cos(it) + cos(angle * pib)) / & +! (cos(solar) +cos(angle * pib)) + edir_t = tile_edir_h * ( cosit + cos(angle) ) / & + ( cossolar + cos(angle) ) + rdir_tmp = edir_t/(tile_edir_h+tile_edif_h) + rtotal_tmp = (edir_t+edif_t)/(tile_edir_h+tile_edif_h) + where( rtotal .le. rth ) + rtotal = rtotal_tmp + rdir = rdir_tmp + end where +! endif !---------------------------------------------------------------- ! calculate normalized BRDF shape function for sloping surface - ann_s = RL_brdf(it_brdf,et_brdf,ra_sl,hb,br,1.0,norm_1,norm_2) + ann_s = RL_brdf(it_brdf,et_brdf,ra_sl,hb,br,1.0,norm_1,norm_2,veclen) + !---------------------------------------------------------------- ! calculate black sky albedo for sloar angle - aa_solars = black_sky(1.0,norm_1,norm_2,it_bk) + aa_solars = black_sky(1.0,norm_1,norm_2,it_bk,veclen) !-------------------------------------------------------------- ! calculate black sky albedo for view angle - aa_views = black_sky(1.0,norm_1,norm_2,et_bk) + aa_views = black_sky(1.0,norm_1,norm_2,et_bk,veclen) !------------------------------------------------------------- aa_slope = & - (rdir * (fv(i, j) * ann_s + (1.0-fv(i, j)) * & - aa_solars) + rdif * (fv(i, j) * aa_views + & - (1.0-fv(i, j)) * aa_white)) / aa_white - - a_eqs = (rtotal - aa_slope) * s_mod(i, j) * & - (1 - s_mod(i, j) * ref_lm(i)) - b_eqs = aa_slope + ref_lm(i) * (1 - aa_slope) * s_mod(i, j) - c_eqs = -ref_lm(i) - - if (abs(a_eqs) .lt. 0.00001) then - ref_bars = -c_eqs / b_eqs - else - ref_bars = (-b_eqs + sqrt(b_eqs**2 -4 * a_eqs * c_eqs)) / & - (2*a_eqs) - endif + (rdir * (tile_fv * ann_s + (1.0-tile_fv) * & + aa_solars) + rdif * (tile_fv * aa_views + & + (1.0-tile_fv) * aa_white)) / aa_white + + a_eqs = (rtotal - aa_slope) * tile_s_mod * & + (1 - tile_s_mod * tile_ref_lm) + b_eqs = aa_slope + tile_ref_lm * (1 - aa_slope) * tile_s_mod +! c_eqs = -ref_lm(i) + +! if (abs(a_eqs) .lt. 0.00001) then +! ref_bars = -c_eqs / b_eqs +! else + ref_bars = (-b_eqs + sqrt(b_eqs**2 +4 * a_eqs * tile_ref_lm)) / & + (2*a_eqs) + where( abs(a_eqs) .lt. 0.00001 ) ref_bars = tile_ref_lm / b_eqs +! endif ref_brdfreals = ann_s * ref_bars / aa_white - ref_terrain(i) = ref_bars * fnn / aa_white - iref_terrain(i, j) = int(ref_terrain(i) * 10000.0 + 0.5) - - if (ref_terrain(i) .ge. 1) then - ref_terrain(i) = 1.0 - iref_terrain(i, j) = int(ref_terrain(i) * 10000.0 + 0.5) - endif + +! ref_terrain(i) = ref_bars * fnn / aa_white +! iref_terrain(i, j) = int(ref_terrain(i) * 10000.0 + 0.5) + tile_ref_terrain = min(ref_bars * fnn / aa_white, 1.0) + tile_iref_terrain = nint(tile_ref_terrain * 10000.0,2) +! if (ref_terrain(i) .ge. 1) then +! ref_terrain(i) = 1.0 +! iref_terrain(i, j) = int(ref_terrain(i) * 10000.0 + 0.5) +! endif ! Should test for these cases in initial tests! (ie must be lt these) ! presently comments as test for ge 90 in initial one ! - if (it_angle(i, j) .ge. 85.0) then - iref_terrain(i, j) = -999 - endif - - if (et_angle(i, j) .ge. 85.0) then - iref_terrain(i, j) = -999 - endif +! if (it_angle(i, j) .ge. 85.0) then +! iref_terrain(i, j) = -999 +! endif +! +! if (et_angle(i, j) .ge. 85.0) then +! iref_terrain(i, j) = -999 +! endif ! if in masked or otherwise unused area you get here ! put in -999 values to indicate null data - else - iref_terrain(i, j) = -999 - endif - else - ref_lm(i) = -999.0 - ref_brdf(i) = -999.0 - ref_terrain(i) = -999.0 - iref_lm(i, j) = -999 - iref_brdf(i, j) = -999 - iref_terrain(i, j) = -999 - endif - enddo +! else +! iref_terrain(i, j) = -999 +! endif + where( it .ge. 85.0*pib .OR. & + et .ge. 85.0*pib .OR. & + tile_mask_self .le. 0 .OR. & + tile_mask_castview .le. 0 .OR. & + tile_mask_castsun .le. 0 ) & + tile_iref_terrain = -999 + else + tile_iref_terrain = -999 + end if + + where (tile_a_mod(1:tilelen) .ge. 0 .and. tile_dn_1(1:tilelen) .gt. 0 ) + iref_lm(istart:iend,j) = tile_iref_lm(1:tilelen) + iref_brdf(istart:iend,j) = tile_iref_brdf(1:tilelen) + iref_terrain(istart:iend,j) = tile_iref_terrain(1:tilelen) + elsewhere +! ref_lm(i) = -999.0 +! ref_brdf(i) = -999.0 +! ref_terrain(i) = -999.0 + iref_lm(istart:iend, j) = -999 + iref_brdf(istart:iend, j) = -999 + iref_terrain(istart:iend, j) = -999 + end where + else + iref_lm(istart:iend, j) = -999 + iref_brdf(istart:iend, j) = -999 + iref_terrain(istart:iend, j) = -999 + end if + enddo enddo -END SUBROUTINE reflectance + END SUBROUTINE reflectance diff --git a/src/Makefile b/src/Makefile index 4fdd2a9..daa97a0 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,7 +1,7 @@ SHELL := /bin/bash -CC := gcc -FC := gfortran -FCFLAGS := -Wall -O3 -fbounds-check +CC := icc +FC := ifort +FCFLAGS := -O3 -check bounds -assume byterecl -xHost -ipo CFLAGS := -O3 EXES := calculate_coefficients \ diff --git a/workflow/nbar.cfg b/workflow/nbar.cfg index e19b62c..ef0e686 100644 --- a/workflow/nbar.cfg +++ b/workflow/nbar.cfg @@ -48,7 +48,7 @@ tc_intermediates = TC_Intermediates rfl_output_dir = Reflectance_Outputs header_slope_target = SLOPE_ANGLE_INPUTS x_tile_size = -1 -y_tile_size = 100 +y_tile_size = 1120 [brdf] factor = 1.0 From e52ac8c6392f684e5638ac4a3e31181ad93629d4 Mon Sep 17 00:00:00 2001 From: Josh Sixsmith Date: Wed, 8 Jun 2016 10:40:21 +1000 Subject: [PATCH 5/5] Convert file to endings to line feed. --- .gitattributes | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitattributes b/.gitattributes index 3b793b2..54110ee 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1,4 +1,4 @@ -* text=auto +* text=lf gaip/version.py export-subst