[GRASS-dev] r.sun and pj_do_proj() calls

GRASS GRASS developer
1

I'm working on translating r.sun to OpenCL to run on the GPU, but I've hit a bit of a snag. I'm hoping y'all can help or at least tell me what I'm looking at.

The root of it is that I can't call pj_do_proj() from within the OpenCL kernel, so I need to pass that in from the outside. For example, on line 1824 of main.c in the current SVN has us calculating the lat/long if we aren't passed in files latin and longin. So for the OpenCL version I'll need to make latin and longin arrays ahead of time instead of calculating it in the loop.

The problem is in rsunlib.c at line 250. How can I avoid this pj_do_proj() call? I was thinking that because I have the lat/long coordinates for all points, I can use that on lines 255 & 256, effectively changing G_projection() to return PROJECTION_LL and avoiding going into the 'if' statement. I don't know if that will work, though. I really don't know what's going on here enough to feel comfortable changing the calculation.

Suggestions? Code?
~Seth

2

Seth Price wrote:

I'm working on translating r.sun to OpenCL to run on the GPU, but I've
hit a bit of a snag. I'm hoping y'all can help or at least tell me
what I'm looking at.

The root of it is that I can't call pj_do_proj() from within the
OpenCL kernel, so I need to pass that in from the outside. For
example, on line 1824 of main.c in the current SVN has us calculating
the lat/long if we aren't passed in files latin and longin. So for the
OpenCL version I'll need to make latin and longin arrays ahead of time
instead of calculating it in the loop.

The problem is in rsunlib.c at line 250. How can I avoid this
pj_do_proj() call? I was thinking that because I have the lat/long
coordinates for all points, I can use that on lines 255 & 256,
effectively changing G_projection() to return PROJECTION_LL and
avoiding going into the 'if' statement. I don't know if that will
work, though. I really don't know what's going on here enough to feel
comfortable changing the calculation.

Suggestions? Code?

I don't know the answer to this. But if no-one else does either, r.sun
should be disabled. An error is preferable to silently producing wrong
answers.

Using different mechanisms for projection in different places looks
really suspicious. If there's a good reason for doing so, it should be
clearly documented. If there isn't a good reason, it should be fixed.
If no-one knows whether it's right or wrong, the module should be
disabled until we do know.

--
Glynn Clements <glynn@gclements.plus.com>

3

On Wed, Jun 9, 2010 at 6:00 PM, Glynn Clements <glynn@gclements.plus.com> wrote:

Seth Price wrote:

I'm working on translating r.sun to OpenCL to run on the GPU, but I've
hit a bit of a snag. I'm hoping y'all can help or at least tell me
what I'm looking at.

The root of it is that I can't call pj_do_proj() from within the
OpenCL kernel, so I need to pass that in from the outside. For
example, on line 1824 of main.c in the current SVN has us calculating
the lat/long if we aren't passed in files latin and longin. So for the
OpenCL version I'll need to make latin and longin arrays ahead of time
instead of calculating it in the loop.

The problem is in rsunlib.c at line 250. How can I avoid this
pj_do_proj() call? I was thinking that because I have the lat/long
coordinates for all points, I can use that on lines 255 & 256,
effectively changing G_projection() to return PROJECTION_LL and
avoiding going into the 'if' statement. I don't know if that will
work, though. I really don't know what's going on here enough to feel
comfortable changing the calculation.

Suggestions? Code?

I don't know the answer to this. But if no-one else does either, r.sun
should be disabled. An error is preferable to silently producing wrong
answers.

For sure. I have added the authors in bcc who may have suggestions.

Using different mechanisms for projection in different places looks
really suspicious. If there's a good reason for doing so, it should be
clearly documented. If there isn't a good reason, it should be fixed.
If no-one knows whether it's right or wrong, the module should be
disabled until we do know.

The first thing to do may be contacting the authors... done with this email.

Markus

4

Does latin/longin work for both locations that pj_do_proj() is called, or just the one in main.c? It looked like some additional calculations were done in the rsunlib.c call.

Would you be able to write some additional code that uses latin/longin for the rsunlib.c call? I fear screwing it up because I'm not sure what it's doing there.

Thanks,
Seth

On Jun 9, 2010, at 9:19 PM, Jaro Hofierka wrote:

Markus Neteler wrote:

On Wed, Jun 9, 2010 at 6:00 PM, Glynn Clements<glynn@gclements.plus.com> wrote:

Seth Price wrote:

I'm working on translating r.sun to OpenCL to run on the GPU, but I've
hit a bit of a snag. I'm hoping y'all can help or at least tell me
what I'm looking at.

The root of it is that I can't call pj_do_proj() from within the
OpenCL kernel, so I need to pass that in from the outside. For
example, on line 1824 of main.c in the current SVN has us calculating
the lat/long if we aren't passed in files latin and longin. So for the
OpenCL version I'll need to make latin and longin arrays ahead of time
instead of calculating it in the loop.

The problem is in rsunlib.c at line 250. How can I avoid this
pj_do_proj() call? I was thinking that because I have the lat/long
coordinates for all points, I can use that on lines 255& 256,
effectively changing G_projection() to return PROJECTION_LL and
avoiding going into the 'if' statement. I don't know if that will
work, though. I really don't know what's going on here enough to feel
comfortable changing the calculation.

Suggestions? Code?

Hi Seth,

Maybe Thomas Huld can answer this better than me, but just a small explanation how it should work.

r.sun effectivelly needs latitude,longitude for every cell. So intention with latin/longin rasters was to get these values in case you don't have a projection defined in grass. If you have a cartesian projection defined in grass, you can get lat/lon on the fly (e.g. using pj_do_proj() or any other method). So if you are able to read the lat/lon values for every cell by other method then you can disable this call.

Jaro

I don't know the answer to this. But if no-one else does either, r.sun
should be disabled. An error is preferable to silently producing wrong
answers.

For sure. I have added the authors in bcc who may have suggestions.

Using different mechanisms for projection in different places looks
really suspicious. If there's a good reason for doing so, it should be
clearly documented. If there isn't a good reason, it should be fixed.
If no-one knows whether it's right or wrong, the module should be
disabled until we do know.

The first thing to do may be contacting the authors... done with this email.

Markus

5

to make things faster (and maybe more consistent), it may just be
better to have lat and long rasters as input to r.sun, instead of
repeatingly calling for a function to convert them.

6

I think that's the best solution for the OpenCL code I'm working on. However, with the normal code, computing the lat/long as needed may take less time, RAM, and hard drive space then pre-computing it. So either situation may be better, depending on how the code is run.
~Seth

On Jun 9, 2010, at 10:16 PM, Yann Chemin wrote:

to make things faster (and maybe more consistent), it may just be
better to have lat and long rasters as input to r.sun, instead of
repeatingly calling for a function to convert them.

7

For reference, this page seems to make a reasonable attempt at explaining the calculations:
http://re.jrc.ec.europa.eu/pvgis/solres/solmod3.htm

~Seth

On Jun 9, 2010, at 9:45 PM, Jaro Hofierka wrote:

Seth Price wrote:

Does latin/longin work for both locations that pj_do_proj() is called,
or just the one in main.c? It looked like some additional calculations
were done in the rsunlib.c call.

In my opinion pj_do_proj() should be called only once, in main.c. I don't see a real reason why it is also in rsunlib.c, so I believe it should be commented out here. If I recall correctly rsunlib.c was created by Thomas so I hope he can confirm this.

Jaro

Would you be able to write some additional code that uses latin/longin
for the rsunlib.c call? I fear screwing it up because I'm not sure what
it's doing there.

Thanks,
Seth

On Jun 9, 2010, at 9:19 PM, Jaro Hofierka wrote:

Markus Neteler wrote:

On Wed, Jun 9, 2010 at 6:00 PM, Glynn
Clements<glynn@gclements.plus.com> wrote:

Seth Price wrote:

I'm working on translating r.sun to OpenCL to run on the GPU, but I've
hit a bit of a snag. I'm hoping y'all can help or at least tell me
what I'm looking at.

The root of it is that I can't call pj_do_proj() from within the
OpenCL kernel, so I need to pass that in from the outside. For
example, on line 1824 of main.c in the current SVN has us calculating
the lat/long if we aren't passed in files latin and longin. So for the
OpenCL version I'll need to make latin and longin arrays ahead of time
instead of calculating it in the loop.

The problem is in rsunlib.c at line 250. How can I avoid this
pj_do_proj() call? I was thinking that because I have the lat/long
coordinates for all points, I can use that on lines 255& 256,
effectively changing G_projection() to return PROJECTION_LL and
avoiding going into the 'if' statement. I don't know if that will
work, though. I really don't know what's going on here enough to feel
comfortable changing the calculation.

Suggestions? Code?

Hi Seth,

Maybe Thomas Huld can answer this better than me, but just a small
explanation how it should work.

r.sun effectivelly needs latitude,longitude for every cell. So
intention with latin/longin rasters was to get these values in case
you don't have a projection defined in grass. If you have a cartesian
projection defined in grass, you can get lat/lon on the fly (e.g.
using pj_do_proj() or any other method). So if you are able to read
the lat/lon values for every cell by other method then you can disable
this call.

Jaro

I don't know the answer to this. But if no-one else does either, r.sun
should be disabled. An error is preferable to silently producing wrong
answers.

For sure. I have added the authors in bcc who may have suggestions.

Using different mechanisms for projection in different places looks
really suspicious. If there's a good reason for doing so, it should be
clearly documented. If there isn't a good reason, it should be fixed.
If no-one knows whether it's right or wrong, the module should be
disabled until we do know.

The first thing to do may be contacting the authors... done with this
email.

Markus

8

I'm looking for when that pj_do_proj() call was added (and by whom and why). I'm looking for the full history of this file:
http://newtrac.osgeo.org/grass/browser/grass/trunk/raster/r.sun/rsunlib.c

I can't seem to go back far enough, though. This is the furthest I can get:
http://newtrac.osgeo.org/grass/browser/grass/trunk/raster/r.sun2/rsunlib.c?rev=34556

Which came from here:
http://newtrac.osgeo.org/grass/log/grass-addons/r.sun_horizon/r.sun2/rsunlib.c?rev=25783

The comments say it came from JRC, but what is that?

I do see that the original function didn't use pj_do_proj(), but I'm not sure why it was added:
http://newtrac.osgeo.org/grass/browser/grass/trunk/raster/r.sun/main.c?rev=38127

Does anyone know what revision pj_do_proj() was added in?
~Seth

On Jun 9, 2010, at 10:16 PM, Yann Chemin wrote:

to make things faster (and maybe more consistent), it may just be
better to have lat and long rasters as input to r.sun, instead of
repeatingly calling for a function to convert them.

9

On Thu, Jun 10, 2010 at 09:10, Seth Price <seth@pricepages.org> wrote:

Which came from here:
http://newtrac.osgeo.org/grass/log/grass-addons/r.sun_horizon/r.sun2/rsunlib.c?rev=25783

The comments say it came from JRC, but what is that?

The JRC seems to be the European Commission Joint Research Centre
http://ec.europa.eu/dgs/jrc/index.cfm

The link to r.sun documentation you gave earlier is a sub-page of that site.

--Wolf

10

On Thu, 10 Jun 2010, Wolf Bergenheim wrote:

On Thu, Jun 10, 2010 at 09:10, Seth Price <seth@pricepages.org> wrote:

Which came from here:
http://newtrac.osgeo.org/grass/log/grass-addons/r.sun_horizon/r.sun2/rsunlib.c?rev=25783

The comments say it came from JRC, but what is that?

The JRC seems to be the European Commission Joint Research Centre
http://ec.europa.eu/dgs/jrc/index.cfm

The link to r.sun documentation you gave earlier is a sub-page of that site.

There is more history still in the old CVS:
http://freegis.org/cgi-bin/viewcvs.cgi/grass/src/raster/r.sun/

I had a look at the code yesterday and the thing I don't really follow is why the reprojection needs to be done in the com_par() function when it has already been done in calculate(). A few lines before the call to pj_do_proj() in com_par(), the cosine is taken of the latitude, which I don't think would make sense if it wasn't already a latitude (but the call to pj_do_proj() suggests it is a northing which needs to be inverse-projected to latitude).

Paul

11

Paul Kelly wrote:

>> http://newtrac.osgeo.org/grass/log/grass-addons/r.sun_horizon/r.sun2/rsunlib.c?rev=25783
>>
>> The comments say it came from JRC, but what is that?
>
> The JRC seems to be the European Commission Joint Research Centre
> http://ec.europa.eu/dgs/jrc/index.cfm
>
> The link to r.sun documentation you gave earlier is a sub-page of that site.

There is more history still in the old CVS:
http://freegis.org/cgi-bin/viewcvs.cgi/grass/src/raster/r.sun/

I had a look at the code yesterday and the thing I don't really follow is
why the reprojection needs to be done in the com_par() function when it
has already been done in calculate(). A few lines before the call to
pj_do_proj() in com_par(), the cosine is taken of the latitude, which I
don't think would make sense if it wasn't already a latitude (but the call
to pj_do_proj() suggests it is a northing which needs to be
inverse-projected to latitude).

The pj_do_proj() in com_par() is an inverse projection, i.e. it
projects from lat/lon back to the location's coordinate system.

This isn't going to be easy to do with the latin/longin arrays. OTOH,
I suspect that it makes the latin/longin arrays pointless; either they
produce identical results to using PROJ, or you get bogus results.

--
Glynn Clements <glynn@gclements.plus.com>

12

Yann Chemin wrote:

> to make things faster (and maybe more consistent), it
> may just be better to have lat and long rasters as input
> to r.sun, instead of repeatingly calling for a function
> to convert them.

it already has that, with the latin= and lonin= options.

It is my intention to remove those two options from r.sun in
trunk, as it turns out that the overhead of reading the
additional two maps and resampling them by libgis is on par
with just doing the proj calc on the fly (at least in my very
few tests). simple_xy location users speak up now if there's
some reason not to do that.

see http://grass.osgeo.org/wiki/R.sun
and https://trac.osgeo.org/grass/ticket/498

(but trac is broken right now due to server migration)

latin/lonin maps needlessly complicate the module + workflow
for very little gain.

Seth:

I think that's the best solution for
the OpenCL code I'm working on. However, with the normal
code, computing the lat/long as needed may take less time,
RAM, and hard drive space then pre-computing it. So either
situation may be better, depending on how the code is run.

If we can figure out a way to do it without those maps it would
be better. :slight_smile: but for now they are there so don't let it be a
blocker for your SoC efforts as the summer is short.

but... I see the input and output projection variables are
swapped in com_par(): the function is *reverse* projecting
from lat/lon back to Cartesian space, presumably so that the
sqrt(a^2+b^2) distance calc is done on axes of equal scale.

but maybe that trig could be fixed up and the dx dy trick to
find the displacement angle improved/replaced..?

(see one of the ps.map wiki page examples where it shows how
to get gedesic convergence angle using proj4 code [implemented
as 'g.region -n' (has memory bug!)] instead of the dx,dy offset
trick)

also note long-standing todo:
"""
Probably the sun position calculation should be replaced with

G_calc_solar_position()

currently used in r.sunmask. G_calc_solar_position() is based on
solpos.c from NREL and very precise.
"""

rsunlib.c is new for "r.sun2"; module history for earlier r.sun
can be found in svn in the grass6.5 raster/r.sun/ dir not the
g65 raster/r.sun2/ dir. (and earlier grass5 history in the old
cvs as pointed out by Paul; in trunk r.sun2/ is renamed r.sun/).

(we are smarter about correct svn merging now, but that one
got past so the change history got split)

regards,
Hamish