Moritz Lennert wrote:
> - A location in EPSG 31370 (Belgian Lambert 1972)
> - Objectif: project the elevation raster layer from the nc_spm_08
> dataset into that Belgian projection dataset
> - A first run of r.proj to get the coordinates of the layer in the new
> projection system. Here's the result:
>
> n=2370907.92051779 s=2378797.24912944 w=-6087464.12326059
> e=-6068211.64408602 rows=1350 cols=1500
>
> As you can see s > n. g.region bails out on that.
>
> But, when you just switch n and s, then everything works fine and the
> elevation dataset gets projected to the correct spot.
>
> Now, I know that this is not the typical reprojection scenario, but I'm
> wondering whether this indicates an underlying issue in r.proj's
> boundary calculation algorithm.
This has nothing to do with r.proj's boundary calculation algorithm,
which isn't used for -p/-g. Instead, it simply projects the south-east
and north-west corners of the input map, and assumes that the
projected points are also the south-east and north-west corners.
Clearly, this will be inaccurate for anything other than projecting
between cylindrical projections and/or lat-lon.
The region-cropping optimisation which is used for projecting (when
the -n isn't given) projects the entire boundary, sampled at the
source resolution, then uses the bounding box of the projected
boundary. AFAICT, this works for all but pathological cases (e.g.
a projection which is tangential to the region being projected).
Can you try the attached patch?
--
Glynn Clements <[email protected]>
Index: raster/r.proj/main.c
===================================================================
--- raster/r.proj/main.c (revision 53270)
+++ raster/r.proj/main.c (working copy)
@@ -354,10 +354,22 @@
G_message(_("Input map <%s@%s> in location <%s>:"),
inmap->answer, setname, inlocation->answer);
+#if 1
+ outcellhd.north = 1e9;
+ outcellhd.south = -1e9;
+ outcellhd.east = 1e9;
+ outcellhd.west = -1e9;
+ bordwalk(&incellhd, &outcellhd, &iproj, &oproj);
+ inorth = outcellhd.north;
+ isouth = outcellhd.south;
+ ieast = outcellhd.east;
+ iwest = outcellhd.west;
+#else
if (pj_do_proj(&iwest, &isouth, &iproj, &oproj) < 0)
G_fatal_error(_("Error in pj_do_proj (projection of input
coordinate pair)"));
if (pj_do_proj(&ieast, &inorth, &iproj, &oproj) < 0)
G_fatal_error(_("Error in pj_do_proj (projection of input
coordinate pair)"));
+#endif
G_format_northing(inorth, north_str, curr_proj);
G_format_northing(isouth, south_str, curr_proj);
_______________________________________________
grass-user mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/grass-user
