[cc'ing the world, hope that's ok]
Seth wrote:
The main problem isn't how
computationally intensive com_par() is, it's that if I want
to rewrite the calculate() main loop into OpenCL, then I
need to rewrite everything that is called by it. When I run
the calculate() OpenCL kernel on the GPU, it can't go back
and call the pj_do_proj() on the host CPU.So we either need to rewrite com_par() to pull pj_do_proj()
outside of the main calculate() loop, or we rewrite all of
pj_do_proj() into OpenCL. The former is difficult, but I
believe the latter is effectively impossible.
I concur.
Thus we need to figure a way for getting pj_do_proj() outside
of the main calculate() loop and passing all it's results into
com_par() in an array. Then we can copy the array to the GPU,
and run the kernel without interruption. Does this make sense?
I've got another idea: all it's doing is setting up a triangle
in Cartesian space to make the trigonometry work easier. I don't
think it is really necessary to do a full transform back to the
map projection to accomplish that, if we neglect the ellipsoid
we can just set east and north axes to nearly the same scale by *cos(lat).
While that might not be as nice as a real map projection, it
actually avoids the 'convergence angle'* problem, which may be
mucking up the right-triangle assumption anyway.
[*] convergence angle is the difference between grid north and
true north at a point, for spearfish this is about 0.9 degrees.
(in grass: `g.region -n`)
try the attached patch against trunk which removes the
pj_do_proj() call and replaces it with a cheap estimation. The
resulting stepsinangle and stepcosangle vars seem within about
0.2-0.4 of the pj_do_proj() version, but after factoring in no
convergence angle error I couldn't tell you how much better or
worse the result is, or vs. the old sunAzimuthAngle method, or
how much this even matters. In my tests I got perfectly identical
output maps before and after the change. Perhaps I was not using
the right combo of options to trigger that, but I was using
shading and no horizon maps, so I'm pretty sure I use that code
fully. ?
bonus: the patch makes r.sun in Mode 2 run 16% faster.
(fwiw I don't really follow why the grid resolution (stepxy) is
in that calc.. shrug)
Honestly, if com_par() was particularly computationally
intensive, that would only be a reason to put as much of it
as possible on the GPU. The CPU is slow.
fwiw my tests were in the spearfish sample dataset:
# DEM (mode 1)
g.region rast=elevation.dem
r.sun -s elevation.dem glob_rad=global_radDay180.14 day=180 time=14
r.colors -e global_radDay180.14 col=grey
d.shadedmap rel=elevation.dem drape=global_radDay180.14 br=30
# Gaussian mound (mode 2)
# r.surf.volcano is from the grass addons repo.
g.region -d
r.surf.volcano out=gauss method=gaussian kurtosis=1
r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 step=0.5
# check for deviations between old and new code:
r.mapcalc "diff = rad.global.ORIG - rad.global.NEW"
r.univar diff
Besides r.sun, my next project will be getting my warping
code integrated into GRASS's warping & interpolation
modules. Depending on how my PhD goes, I may be putting
effort into other GRASS modules this fall, though.What's
related to land cover classification if I'm writing a new
classifier module already?
umm, pass. beyond my area of expertise. [cc list]
Once I get some appropriate 4+ band imagery, I'll be making
some tweaks for the CPU code to use SSE also. Debugging
OpenCL is hell, but it's pretty powerful & versatile
once you get it running. (Looks like I have some 6-band
imagery, I know what I'm doing today...)
From the GRASS download page you can find the Spearfish Imagery
extra mapset, and also the new North Carolina sample dataset.
Both include multiband LANDSAT data. Also the NC dataset has
geotiff versions available for download if you want to push it
straight into gdalwarp. otherwise just use `r.out.gdal` to
create the geotiffs from the grass data.
Also you can go to
http://grass.osgeo.org/wiki/Global_datasets#LANDSAT
where there is a link to where you can access the USGS LANDSAT
archive. (requests may take a few days to process)
regards,
Hamish
(attachments)
rsun_no_pj.diff (2.4 KB)