[GRASS-user] a bug in v.proj?

[Ivan Shmakov]

>>>>> "MN" == Markus Neteler <neteler at itc.it> writes:

[...]

 >>> I seem not to understand the hack (why +ellps=sphere in -t_srs?),
 >>> but it, indeed, works.

 >>  This is related I guess:
 >> http://proj.maptools.org/faq.html#sphere_as_wgs84

 > I don't think that it is related.  I think that it is a GDAL
 > issue. You may check and join my bug report
 > http://trac.osgeo.org/gdal/ticket/1239

[...]

	Well, I've first encountered the ``MODIS L3'' problem when I've
	tried to import several files into GRASS using `r.in.gdal'.
	After a few vector maps were reprojected into the location as
	well, it was noticed, that either rasters or vectors were
	shifted (or both.)

	Now, I think that it's an issue with `v.proj', while `r.in.gdal'
	actually works fine.  Please consider the following (tested on
	GRASS 6.0.2 debian 6.)

	Let me introduce two locations first:

world-ll $ g.proj -p
-PROJ_INFO-------------------------------------------------
name       : Latitude-Longitude
datum      : wgs84
towgs84    : 0.000,0.000,0.000
proj       : ll
ellps      : wgs84
-PROJ_UNITS------------------------------------------------
unit       : degree
units      : degrees
meters     : 1.0
world-ll $ 

modis-sinu $ g.proj -p
-PROJ_INFO-------------------------------------------------
name       : Sinusoidal (Sanson-Flamsteed)
towgs      : 0,0,0
proj       : sinu
R          : 6371007.181
a          : 6371007.181
b          : 6371007.181
es         : 0.0000000000
f          : 0.0000000000
lat_0      : 0.0000000000
lon_0      : 0.0000000000
-PROJ_UNITS------------------------------------------------
unit       : meter
units      : meters
meters     : 1.0
modis-sinu $ 

	... And a `points' file (longitude-latitude coordinates):

world-ll $ cat points-ll.dat
80.00|50.00|a point on L3 (22, 3) tile bottom edge
83.75|53.36|Barnaul
90.00|60.00|a point on L3 (22, 3) tile top edge
world-ll $ v.in.ascii input=points-ll.dat \
               columns='x double precision, y double precision, name varchar(40)' \
               output=marked_points
...
world-ll $ v.out.ascii input=marked_points
80|50|1
83.75|53.36|2
90|60|3
world-ll $ 

	The points on the edges of an L3 tile are to allow the `v.proj'
	result to be easily checked afterwards.  (Rows of the L3 tiles
	are the stripes between consequent 10 degrees parallels, so row
	#3 is bounded by 60dN and 50dN.)  The ``Barnaul'' point is to be
	checked visually.

	Then, I reproject the map:

modis-sinu $ v.proj input=marked_points location=world-ll mapset=ivan
Input Projection Parameters: +proj=latlong +a=6378137 +rf=298.257223563 +no_defs +towgs84=0.000,0.000,0.000
Input Unit Factor: 1

Output Projection Parameters: +proj=sinu +R=6371007.181 +lat_0=0.0000000000 +lon_0=0.0000000000 +no_defs
Output Unit Factor: 1
...
modis-sinu $ v.out.ascii input=marked_points
5740503.95315973|5538668.85235319|1
5581642.24347567|5912856.23809764|2
5029005.65975492|6653142.01247702|3
modis-sinu $ 

	The result seems not to be correct.  The L3 tiles in sinusoidal
	projection are the squares comprised of 1200 x 1200 pixels each,
	and each pixel is a 926.62543306 meters square.  There are 36 x
	18 L3 tiles covering the whole Earth surface.

	Rows of the tiles are numbered north-to-south, with 90..80dN row
	numbered 0th.  So, 3rd row is 60..50dN, and it corresponds to
	Y-range of:

(- (* 926.62543306 1200 9) (* 926.62543306 1200 3))
;; 60dN
=> 6671703.118031999
(- (* 926.62543306 1200 9) (* 926.62543306 1200 4))
;; 50dN
=> 5559752.598359999

	While the `v.proj' results in the points being shifted by the
	amount of 19..21 km:

(- 6653142.01247702 6671703.118031999)
;; 60dN
=> -18561.105554979295
(- 5538668.85235319 5559752.598359999)
;; 50dN
=> -21083.746006809175

	When I use `proj' to ``reproject'' the data, it seems to work
	just fine:

modis-sinu $ sed -e 's,^\([^|]*\)|\([^|]*\)\(|.*\),\1 \2 \3,' \
                 < points-ll.dat \
                 | proj $(g.proj -jf) \
                 | sed -e 's,^\([^[:blank:]]*\)[[:blank:]]\+\([^[:blank:]]*\)[[:blank:]]\+\(|.*\),\1|\2\3,' \
                 | tee points-sinu.dat
5717984.13|5559752.60|a point on L3 (22, 3) tile bottom edge
5557613.28|5933367.97|Barnaul
5003777.34|6671703.12|a point on L3 (22, 3) tile top edge
modis-sinu $ v.in.ascii input=points-sinu.dat \
                 columns='x double precision, y double precision, name varchar(40)' \
                 output=marked_points_proj
...
modis-sinu $ 

	The points are visually at the top and at the bottom of the
	imported raster, and ``Barnaul'' seems to be at the right place.