r.sun segfault in mode 2

[neteler]

With the NC data set, this crashes:

GRASS 6.4.2svn (nc_spm_08):~ > g.region rast=elevation
GRASS 6.4.2svn (nc_spm_08):~ > r.sun elevin=elevation aspin=aspect slopein=slope lat=36.107 day=75 beam_rad=b_rad.075 diff_rad=d_rad.075 refl_rad=r_rad.075 
Mode 2: integrated daily irradiation for a given day of the year
Segmentation fault

Debugging output which looks like bad initialization:

GRASS 6.4.2svn (nc_spm_08):~ > gdb r.sun
GNU gdb (GDB) 7.1-1mdv2010.1 (Mandriva Linux release 2010.1)
...
(gdb) r elevin=elevation aspin=aspect slopein=slope lat=36.107 day=75 beam_rad=b_rad.075 diff_rad=d_rad.075 refl_rad=r_rad.075
Starting program: /home/neteler/grass64/dist.x86_64-unknown-linux-gnu/bin/r.sun elevin=elevation aspin=aspect slopein=slope lat=36.107 day=75 beam_rad=b_rad.075 diff_rad=d_rad.075 refl_rad=r_rad.075

[Thread debugging using libthread_db enabled]
Mode 2: integrated daily irradiation for a given day of the year
   0%
Program received signal SIGSEGV, Segmentation fault.
pj_transform (srcdefn=0x0, dstdefn=0x0, point_count=1, point_offset=0, x=0x7fffffffd0e0, y=0x7fffffffd0d8, 
    z=0x7fffffffd0d0) at pj_transform.c:98
98          if( srcdefn->is_geocent )

(gdb) bt full
#0  pj_transform (srcdefn=0x0, dstdefn=0x0, point_count=1, point_offset=0, x=0x7fffffffd0e0, y=0x7fffffffd0d8, 
    z=0x7fffffffd0d0) at pj_transform.c:98
        i = <value optimized out>
#1  0x00007ffff7bd9b26 in pj_do_proj (x=0x7fffffffd268, y=0x7fffffffd270, info_in=0x60d760, info_out=0x60d7e0)
    at do_proj.c:100
        ok = 1092827616
        u = 0
        v = 0
        h = 0
#2  0x00000000004080e1 in calculate (singleSlope=1.2135880960744616e-311, singleAspect=3.6562339989189768e-319, 
    singleAlbedo=0.20000000000000001, singleLinke=3, gridGeom=...) at main.c:1865
        i = 1
        j = 1
        l = 1350
        someRadiation = 1
        numRows = 1350
        arrayOffset = 1
        lum = 6.9533488366228911e-310
        q1 = 3.1829936871890069e-311
        dayRad = 3.1514668911888089e-317
        latid_l = 6.9533558072950894e-310
        cos_u = 2.0779926761062822e-317
        cos_v = 6.9533558073061565e-310
        sin_u = 6.953349070413816e-310
        sin_v = 0
        sin_phi_l = 0
        tan_lam_l = 0
        zmax = 156.32986450195312
        longitTime = 0
        locTimeOffset = 1.4821969375237396e-323
        latitude = 215010
        longitude = 630010
        coslat = 630000
        sunGeom = {lum_C11 = 2.0969134140794186e-317, lum_C13 = 6.3691388397870058e-314, 
          lum_C22 = 6.9533488366228911e-310, lum_C31 = 6.9533491751815725e-310, lum_C33 = 6.9533558072603072e-310, 
          sunrise_time = 2.1219957904712067e-314, sunset_time = 6.9533558072603072e-310, 
          timeAngle = 4.9406564584124654e-324, sindecl = 0.030294190998040194, cosdecl = 0.99954102566716807}
        sunVarGeom = {isShadow = -11720, z_orig = 127.38738250732422, zmax = 156.32986450195312, zp = 127.38738250732422, 
          solarAltitude = 0, sinSolarAltitude = 0, tanSolarAltitude = 6.9533488354035371e-310, 
          solarAzimuth = 6.9533490691820115e-310, sunAzimuthAngle = 6.9531436082559572e-310, stepsinangle = 0, 
          stepcosangle = 0}
        sunSlopeGeom = {longit_l = 6.9533558072535879e-310, lum_C31_l = 2.1219957909652723e-314, 
          lum_C33_l = 6.9533488359529381e-310, slope = 0.075450657666172191, aspect = 5.0632796287624782}
        sunRadVar = {cbh = 1, cdh = 1, linke = 3, G_norm_extra = 1381.7904736323146, alb = 0.20000000000000001}
#3  0x000000000040488b in main (argc=9, argv=0x7fffffffd638) at main.c:768
        singleSlope = 1.2135880960744616e-311
        singleAspect = 3.6562339989189768e-319
        singleAlbedo = 0.20000000000000001
        singleLinke = 3
        module = 0x7ffff7bd4f40
        parm = {elevin = 0x7ffff7bd4ea0, aspin = 0x613da0, aspect = 0x613e70, slopein = 0x613f40, slope = 0x614010, 
          linkein = 0x6140e0, lin = 0x6141b0, albedo = 0x614280, longin = 0x6145c0, alb = 0x614350, latin = 0x614420, 
          lat = 0x6144f0, coefbh = 0x614690, coefdh = 0x614760, incidout = 0x6149d0, beam_rad = 0x614aa0, 
          insol_time = 0x614b70, diff_rad = 0x614c40, refl_rad = 0x614d10, glob_rad = 0x614de0, day = 0x614eb0, 
---Type <return> to continue, or q <return> to quit---
          step = 0x614f80, declin = 0x615050, ltime = 0x615120, dist = 0x6151f0, horizon = 0x614830, 
          horizonstep = 0x614900, numPartitions = 0x6152c0, civilTime = 0x615390}
        flag = {shade = 0x7ffff7bd4e60, saveMemory = 0x615490}
        gridGeom = {xp = 6.9531436082559572e-310, yp = 0, xx0 = 6.9533491762080433e-310, yy0 = 0, 
          xg0 = 6.9533491624931762e-310, yg0 = 6.9533491654678467e-310, stepx = 10, stepy = 10, deltx = 15000, 
          delty = 13500, stepxy = 10, sinlat = 6.9533558072840223e-310, coslat = 6.9533558072852081e-310}

[hamish]

Hi,

I can reproduce with the 6.4.0~rc6+42329-3 (essentially 6.4.0) Debian/stable package on amd64.

the problem is due to the lat= option.

note in grass65+trunk I've removed it (2 years ago) due to incompatibility with rsunlib, but in practice it was not much help anyway.

GRASS65> r.sun --help
...
     lat   This does nothing. It is retained for backwards compatibility

see ​https://trac.osgeo.org/grass/ticket/498#comment:14

I think the bug here is that I never backported the fix to the relbr64.

Hamish

ps- see also in #498 comment 8+ (and attachment image) and maybe the r.sun wiki page about the difference that slopein= and aspin= maps make, in fact they ruin the result! until that is fixed, you get better output if you just leave those two options out and let the slope/aspect calculations happen as the module runs.

[hamish]

r.sun2 in 6.4svn sync'd with existing 6.5svn fixes in r47993.

No more segfault, closing ticket.

Hamish