Folks,
I committed changes to datum_shift.c to test for null values of the datum names. I also committed changes to correct the "dubious behavior" that Eric pointed out in r.proj and did the same for v.proj and s.proj.
I compiled the code and ran the updated r.proj without problems. I have not run v.proj or s.proj. Further, all my test locations have their datum and ellipse specified, so they don't actually test the new code. In brief, more testing is advisable.
Roger Miller
On Fri, May 24, 2002 at 07:09:29PM -0400, rgrmill wrote:
Folks,
I committed changes to datum_shift.c to test for null values of the datum
names. I also committed changes to correct the "dubious behavior" that
Eric pointed out in r.proj and did the same for v.proj and s.proj.I compiled the code and ran the updated r.proj without problems. I have
not run v.proj or s.proj. Further, all my test locations have their datum
and ellipse specified, so they don't actually test the new code. In
brief, more testing is advisable.
Roger,
I have tried r.proj again. It still crashes for a standard GRASS lat/long
location (when generating on GRASS startup).
Contents are:
cat grassdata/europa/PERMANENT/PROJ_INFO
name: Latitude-Longitude
proj: ll
ellps: wgs84
When adding a map datum with g.setproj:
cat grassdata/europa/PERMANENT/PROJ_INFO
name: Latitude-Longitude
datum: wgs84
dx: 0.000000
dy: 0.000000
dz: 0.000000
proj: ll
ellps: wgs84
r.proj no longer crashes. For backward compatibility I suggest to also
support the case when the map datum is absent.
Thanks,
Markus
On Fri, May 24, 2002 at 07:09:29PM -0400, rgrmill wrote:
Folks,
I committed changes to datum_shift.c to test for null values of the datum
names. I also committed changes to correct the "dubious behavior" that
Eric pointed out in r.proj and did the same for v.proj and s.proj.I compiled the code and ran the updated r.proj without problems. I have
not run v.proj or s.proj. Further, all my test locations have their datum
and ellipse specified, so they don't actually test the new code. In
brief, more testing is advisable.
Ah - it seems that the update was only commited to the release branch
but not HEAD. Please
- first commit to HEAD,
- the commit to release branch, if appropriate.
Thanks,
Markus
On Mon, May 27, 2002 at 11:00:30AM +0200, Markus Neteler wrote:
On Fri, May 24, 2002 at 07:09:29PM -0400, rgrmill wrote:
>
> Folks,
>
> I committed changes to datum_shift.c to test for null values of the datum
> names. I also committed changes to correct the "dubious behavior" that
> Eric pointed out in r.proj and did the same for v.proj and s.proj.
>
> I compiled the code and ran the updated r.proj without problems. I have
> not run v.proj or s.proj. Further, all my test locations have their datum
> and ellipse specified, so they don't actually test the new code. In
> brief, more testing is advisable.Ah - it seems that the update was only commited to the release branch
but not HEAD. Please
- first commit to HEAD,
- the commit to release branch, if appropriate.
... while
Author: roger
Update of /grassrepository/grass/src/libes/proj
In directory doto:/tmp/cvs-serv27797
Modified Files:
datum_shift.c
Log Message:
Added checks for null or empty datum names
was submitted to HEAD and not to release branch as well.
Roger, please verify for future submissions. I'll fix that now
in CVS.
Markus
On Fri, May 24, 2002 at 07:09:29PM -0400, rgrmill wrote:
Folks,
I committed changes to datum_shift.c to test for null values of the datum
names. I also committed changes to correct the "dubious behavior" that
Eric pointed out in r.proj and did the same for v.proj and s.proj.I compiled the code and ran the updated r.proj without problems. I have
not run v.proj or s.proj. Further, all my test locations have their datum
and ellipse specified, so they don't actually test the new code. In
brief, more testing is advisable.
Roger,
the v.proj seems to work now without datum/ellipse (tested for lat/long)
while r.proj continues to crash.
Please test the code again, you could disable datum/ellipse in
PROJ_INFO with # easily.
Markus
The other way around, please:
If critical: commit to release branch
then commit to HEAD branch
(make sure it works!)
If not critical: commit to HEAD branch
On Mon, May 27, 2002 at 11:00:30AM +0200, Markus Neteler wrote:
On Fri, May 24, 2002 at 07:09:29PM -0400, rgrmill wrote:
Ah - it seems that the update was only commited to the release branch
but not HEAD. Please
- first commit to HEAD,
- the commit to release branch, if appropriate.
On Mon, May 27, 2002 at 04:07:59PM +0200, Bernhard Reiter wrote:
The other way around, please:
If critical: commit to release branch
then commit to HEAD branch
(make sure it works!)
Well, since this was a new feature and not yet fully functional,
I recommend to "play" in the HEAD. At time release branch
is broken... (HEAD as well). But I think that only fully functional
uploads to release branch make sense.
Second problem is that part was commited to HEAD, other part
to release branch.
Markus
On Mon, May 27, 2002 at 04:11:44PM +0200, Markus Neteler wrote:
On Mon, May 27, 2002 at 04:07:59PM +0200, Bernhard Reiter wrote:
> The other way around, please:
>
> If critical: commit to release branch
> then commit to HEAD branch
> (make sure it works!)Well, since this was a new feature
Then it should not get uploaded on the release branch at all!
Only release critical fixes, please.
and not yet fully functional,
I recommend to "play" in the HEAD. At time release branch
is broken... (HEAD as well). But I think that only fully functional
uploads to release branch make sense.
True. But avoid breaking the release branch if possible.
Second problem is that part was commited to HEAD, other part
to release branch.
I agree that this should not happen.
Good that you've noticed.
... sorry to insist... Is there any hope to get
[rvs].proj functional again and to release 5.0.0pre5?
I am not sure if I should try to fix Roger's fixes 
(maybe not).
Roger, if you commented the 'datum' entry in PROJ_INFO
(both, and/or target and source location), you could test
if it works for you (to simulate a non-US location).
Thanks,
Markus
On Tue, May 28, 2002 at 04:58:06PM +0200, Markus Neteler wrote:
... sorry to insist... Is there any hope to get
[rvs].proj functional again and to release 5.0.0pre5?
We have to.
Can you write another short summary on the current state so that maybe
another developer can help out here?
On Tue, May 28, 2002 at 04:58:06PM +0200, Markus Neteler wrote:
... sorry to insist... Is there any hope to get
[rvs].proj functional again and to release 5.0.0pre5?I am not sure if I should try to fix Roger's fixes
(maybe not).
Roger, if you commented the 'datum' entry in PROJ_INFO
(both, and/or target and source location), you could test
if it works for you (to simulate a non-US location).
Think I found the bug in r.proj -- a type typo. The same error didn't
appear in v.proj or s.proj.
I'm having problems with UTM to Albers here. With r.proj, I get a
message about shifting datum NAD83->NAD27 (this should be correct), but
then it tells me the Input Map is outside the current region (this is
incorrect). The error occurs in bordwalk where it checks the
boundaries. There's a comparison between some xmin, ymin, etc. and the
region boundaries, and it comes up that ymin is slightly larger (1E-6)
than to_hd->north (this difference is way way below the significant
decimal places of the projection, which isn't better than a couple
places). There's some subtle numerical bug in there...
--
Eric G. Miller <egm2@jps.net>
Eric G. Miller wrote:
I'm having problems with UTM to Albers here. With r.proj, I get a
message about shifting datum NAD83->NAD27 (this should be correct), but
then it tells me the Input Map is outside the current region (this is
incorrect). The error occurs in bordwalk where it checks the
boundaries. There's a comparison between some xmin, ymin, etc. and the
region boundaries, and it comes up that ymin is slightly larger (1E-6)
than to_hd->north (this difference is way way below the significant
decimal places of the projection, which isn't better than a couple
places). There's some subtle numerical bug in there...
If "ymin is slightly larger (1E-6) than to_hd->north" then, unless
these values have been miscalculated, that implies that the source and
destination regions almost intersect, but not quite.
IOW, if you moved the border slightly, you would get some non-null
cells in the output but, as things stand, the output will consist
entirely of null values.
[Aside: there are cases where the "border walking" algorithm won't
produce the right results (e.g. converting polar stereographic to any
cylindrical projection), but I doubt that such cases are relevant to
"real world" use.]
--
Glynn Clements <glynn.clements@virgin.net>
On Wed, 29 May 2002, Glynn Clements wrote:
> boundaries. There's a comparison between some xmin, ymin, etc. and the
> region boundaries, and it comes up that ymin is slightly larger (1E-6)
> than to_hd->north (this difference is way way below the significant
> decimal places of the projection, which isn't better than a couple
> places). There's some subtle numerical bug in there...
I don't think this has anything to do with the bug you are looking for,
but:
No, not a numerical bug. At the start of the bordwalk routine xmin, xmax,
ymin, ymax are intentionally set to (unreasonable) values slightly
larger/smaller than east, west, north, south respectively. At the end, if
any of these variables are still outside the borders it means no point
along borders of source map fell inside the destination area.
So, if you get ymin > to_hd->north at the end of bordwalk, it means all of
the cell centers along the lower border of the source map had (projected)
y-values greater than 'north' in the destination area. Thus 'input area is
outside current region'.
[Aside: there are cases where the "border walking" algorithm won't
produce the right results (e.g. converting polar stereographic to any
cylindrical projection), but I doubt that such cases are relevant to
"real world" use.]
Have you tested this in practice? With bordwalk r.proj should be able to
produce maps composed of discontinous areas. I tested it with some other
projections when I wrote the code, but not with polar stereographic. One
reason for writing bordwalk.c was that in the previous version of r.proj
you could only do reliably cylindrical->cylindrical projections.
Morten Hulden
On Wed, 29 May 2002, Morten Hulden wrote:
So, if you get ymin > to_hd->north at the end of bordwalk, it means all of
the cell centers along the lower border of the source map had (projected)
y-values greater than 'north' in the destination area. Thus 'input area is
outside current region'.
Small correction: actually it means _none_ of the border cells, regardless
of the edge, had projected y-values smaller than to_hd->north.
> [Aside: there are cases where the "border walking" algorithm won't
> produce the right results (e.g. converting polar stereographic to any
> cylindrical projection), but I doubt that such cases are relevant to
> "real world" use.]Have you tested this in practice? With bordwalk r.proj should be able to
produce maps composed of discontinous areas. I tested it with some other
projections when I wrote the code, but not with polar stereographic. One
reason for writing bordwalk.c was that in the previous version of r.proj
you could only do reliably cylindrical->cylindrical projections.
Addition: proj with fail with an error if you try to project a point that
is a reference point for the projection in question (resulting in division
by zero usually). I believe that in polar projections 0,90 and 0,-90 are
such points. I also believe [rvs].proj in such cases will abort. Maybe
this should be changed to just ignoring the point/cell that could not be
projected instead of aborting?
Morten Hulden
Morten Hulden wrote:
> [Aside: there are cases where the "border walking" algorithm won't
> produce the right results (e.g. converting polar stereographic to any
> cylindrical projection), but I doubt that such cases are relevant to
> "real world" use.]Have you tested this in practice? With bordwalk r.proj should be able to
produce maps composed of discontinous areas. I tested it with some other
projections when I wrote the code, but not with polar stereographic. One
reason for writing bordwalk.c was that in the previous version of r.proj
you could only do reliably cylindrical->cylindrical projections.
If I disable the "special cases" code, the obvious UPS->lat/lon test
fails in the expected manner (i.e. I don't get anything north of the
circle inscribing the region). IOW, the actual "border walking" part
of bordwalk() doesn't suffice in this case.
As the "special cases" code only tests four specific points, it won't
always help; although it will help when projecting UPS to the entire
lat/lon range, as one of those four points happens to be the north
pole (minus half a cell).
However, as I said previously, I doubt that these cases are of
real-world importance; such transformations inherently have extreme
distortion, so the projected raster isn't of much use even when it is
complete.
--
Glynn Clements <glynn.clements@virgin.net>
Morten Hulden wrote:
> > [Aside: there are cases where the "border walking" algorithm won't
> > produce the right results (e.g. converting polar stereographic to any
> > cylindrical projection), but I doubt that such cases are relevant to
> > "real world" use.]
>
> Have you tested this in practice? With bordwalk r.proj should be able to
> produce maps composed of discontinous areas. I tested it with some other
> projections when I wrote the code, but not with polar stereographic. One
> reason for writing bordwalk.c was that in the previous version of r.proj
> you could only do reliably cylindrical->cylindrical projections.Addition: proj with fail with an error if you try to project a point that
is a reference point for the projection in question (resulting in division
by zero usually). I believe that in polar projections 0,90 and 0,-90 are
such points. I also believe [rvs].proj in such cases will abort. Maybe
this should be changed to just ignoring the point/cell that could not be
projected instead of aborting?
When computing the region of interest, cells which fail to project
should just be ignored. When performing the actual projection, cells
which fail to project should be NULL.
Also, I've remembered another issue with r.proj:
When projecting lat/lon to pseudo-cylindrical projections where the
entire globe fits within the region, the resulting map tends to
"wrap"; i.e. the area outside of the "ellipse" contains non-null
cells.
--
Glynn Clements <glynn.clements@virgin.net>
On Wed, 29 May 2002, Glynn Clements wrote:
When computing the region of interest, cells which fail to project
should just be ignored. When performing the actual projection, cells
which fail to project should be NULL.
Should be fairly easy to fix: Find the places in the code where
[src].proj exits with 'Error in do_proj' and change action according to
above.
(I don't have a GRASS installation to work with right now, but maybe
someone else is willing to take on this task)
Also, I've remembered another issue with r.proj:
When projecting lat/lon to pseudo-cylindrical projections where the
entire globe fits within the region, the resulting map tends to
"wrap"; i.e. the area outside of the "ellipse" contains non-null
cells.
Can you confirm this happens with recent versions of proj
(post-December-2001).
The way r.proj works (omitting here the initial testing/trimming
done in bordwalk.c) is:
For each cell in the destination region
-project the cell center into source map
-find the category value of the nearest cell
(alternatively using cubic etc interpolation)
-assign that value to the cell in the destination region
For a pseudo-cylindrical destination there will a number of cells outside
the ellipse that should become NULL. (The real 'borders' in this case is
actually the ellipse, not the window borders given in DEFAULT_WIND).
Projecting a co-ordinate outside the ellipse should result in proj
returning * (inf?) for either lon or lat or both.
The fix would be to look up the place in the code where r.proj checks the
return values from proj and if either lat or lon is invalid, assign the
destination cell a NULL value without even trying to find the category
value of the nearest cell in the source map.
I faintly remember r.proj already does something like this, although the
checks may be incomplete (maybe it's only checking for NULL cells in the
source map), but not having the code right now it is impossible for me to
say.
Morten Hulden
Morten Hulden wrote:
> When computing the region of interest, cells which fail to project
> should just be ignored. When performing the actual projection, cells
> which fail to project should be NULL.Should be fairly easy to fix: Find the places in the code where
[src].proj exits with 'Error in do_proj' and change action according to
above.(I don't have a GRASS installation to work with right now, but maybe
someone else is willing to take on this task)
OK.
> Also, I've remembered another issue with r.proj:
>
> When projecting lat/lon to pseudo-cylindrical projections where the
> entire globe fits within the region, the resulting map tends to
> "wrap"; i.e. the area outside of the "ellipse" contains non-null
> cells.Can you confirm this happens with recent versions of proj
(post-December-2001).
This is with the current GRASS CVS head.
The way r.proj works (omitting here the initial testing/trimming
done in bordwalk.c) is:For each cell in the destination region
-project the cell center into source map
-find the category value of the nearest cell
(alternatively using cubic etc interpolation)
-assign that value to the cell in the destination regionFor a pseudo-cylindrical destination there will a number of cells outside
the ellipse that should become NULL. (The real 'borders' in this case is
actually the ellipse, not the window borders given in DEFAULT_WIND).
Projecting a co-ordinate outside the ellipse should result in proj
returning * (inf?) for either lon or lat or both.
That clearly isn't the case. If it were, the cells outside of the
ellipse would be garbage, but they are actually "repeats" of the data
inside the ellipse (albeit more distorted), which seems to suggest
periodicity.
In this test case (Wagner I), the repeat happens in both latitude and
longitude, but the latitudinal repeat is a mirror image in both
directions (i.e. 180° rotation).
Anyway, I'll look into it further.
--
Glynn Clements <glynn.clements@virgin.net>
On Wed, May 29, 2002 at 12:27:23AM -0700, Eric G. Miller wrote:
I'm having problems with UTM to Albers here. With r.proj, I get a
Never mind. Something is horked with a data layer I was using to set
the output region extents...
--
Eric G. Miller <egm2@jps.net>
Glynn Clements wrote:
> > When computing the region of interest, cells which fail to project
> > should just be ignored. When performing the actual projection, cells
> > which fail to project should be NULL.
>
> Should be fairly easy to fix: Find the places in the code where
> [src].proj exits with 'Error in do_proj' and change action according to
> above.
>
> (I don't have a GRASS installation to work with right now, but maybe
> someone else is willing to take on this task)OK.
OK, done.
> The way r.proj works (omitting here the initial testing/trimming
> done in bordwalk.c) is:
>
> For each cell in the destination region
> -project the cell center into source map
> -find the category value of the nearest cell
> (alternatively using cubic etc interpolation)
> -assign that value to the cell in the destination region
>
> For a pseudo-cylindrical destination there will a number of cells outside
> the ellipse that should become NULL. (The real 'borders' in this case is
> actually the ellipse, not the window borders given in DEFAULT_WIND).
> Projecting a co-ordinate outside the ellipse should result in proj
> returning * (inf?) for either lon or lat or both.That clearly isn't the case. If it were, the cells outside of the
ellipse would be garbage, but they are actually "repeats" of the data
inside the ellipse (albeit more distorted), which seems to suggest
periodicity.In this test case (Wagner I), the repeat happens in both latitude and
longitude, but the latitudinal repeat is a mirror image in both
directions (i.e. 180° rotation).Anyway, I'll look into it further.
Done. The following patch produces the "desired" result (i.e. the way
that pseudo-cylindrical world maps normally look in illustrations),
but is it the right thing to do?
diff -u -r1.2 PJ_urmfps.c
--- src/libes/proj/PJ_urmfps.c 20 Apr 2002 19:13:44 -0000 1.2
+++ src/libes/proj/PJ_urmfps.c 30 May 2002 04:30:08 -0000
@@ -17,8 +17,12 @@
}
INVERSE(s_inverse); /* sphere */
xy.y /= P->C_y;
+ if (fabs(xy.y) > HALFPI)
+ I_ERROR
lp.phi = aasin(sin(xy.y) / P->n);
lp.lam = xy.x / (C_x * cos(xy.y));
+ if (fabs(lp.lam) > PI)
+ I_ERROR
return (lp);
}
FREEUP; if (P) pj_dalloc(P); }
If it is, instinct suggests that the same issue will apply to many
(most?) of the projections, not just Wagner I.
Also, I had to disable the remaining call to bordwalk() in r.proj,
otherwise it cropped the result. Realistically, I think that we need a
flag to allow the bordwalk() call(s) to be disabled for the cases
where they don't work correctly.
In the longer term, we should re-think the architecture of the PROJ
library. There are many practical applications where it would be
useful to have more information about the overall "nature" of a
projection than just a pair of functions which project individual
points.
--
Glynn Clements <glynn.clements@virgin.net>