[GRASS-user] Tiling gtopo30

[Zenon Panoussis]

Hello everyone

I am trying to tile the global gtopo30 data into somewhat manageable
rasters. After a lot of googling and trial and error, I'm still not
further than error. This is what I've done:

                                     90
I created a new location, wgs84, -180  180,  00:00:30, 00:00:30,
                                    -60

then run

g.region -d
r.in.bin -bs input=/vol/gis/gtopo30/W180N90.DEM output=gtopo30_W180N90 \
         bytes=2 north=90 south=40 east=-140 west=-180 r=6000 c=4800
# and so on for all 27 non-Antarctic parts

r.patch in=gtopo30_E020N40,gtopo30_E020N90,gtopo30_E020S10[,...] # all 27 parts
g.region rast=gtopo30_full_noant
r.null gtopo30_full_noant setnull=-9999

So far, so well. 'd.rast gtopo30_full_noant' shows the world except
Antarctis with blue lowlands, purple mountains and white oceans.
It's after this point that things go wrong.

v.mkgrid position=region map=gtopo30_grid grid=30,72

for i in `seq 1 2160`
do
  while [ ${#i} -lt 4 ]
    do i=0$i
  done

  g.region rast=gtopo30_full_noant
  v.extract in=gtopo30_grid out=gtopo30_grid_$i list=$i
  g.region vect=gtopo30_grid_$i
  r.mapcalc "gtopo30_$i = gtopo30_full_noant"
#  g.region rast=gtopo30_$i (makes no difference, is same as the vector)
  r.out.tiff -t -v in=gtopo30_$i out=gtopo30_$i
  g.remove vect=gtopo30_grid_$i
  g.remove rast=gtopo30_$i
done

This produces 2160 tile tiffs of which 90 have a size (1-5 MB) and
the rest are 168 bytes each, obviously nothing but a header. That's
clearly wrong, land covers more than 4% of the planet. I finish off
anyway with 'gdaltindex gtopo30_idx.shp /path/to/gtopo30_*.tif' and
feed the resulting .shp to the UMN mapserver. That gives me a stripe
of colour approximately 30 degrees wide covering east Siberia, the
Koreas, the Philippines and part of Australia. If I continue east,
I find the entire world between 180E and 380E, minus the Korea stripe.

It looks like a region error, but I can't see a region error in the
script. Can anyone else? Perhaps I should mention that mapserver has
EXTENT -180 -90 180 90 and ETOPO2 and all my other layers show perfectly
where they should, so the error is not in mapserver.

Another problem is that I see bands of different elevation where the
elevation is the same. So, for instance, the north of the Caspian sea
shows at 0-99m and the south of it at 100-199m. The bands are about
5 degrees of latitude wide all the way from east to west and they are
repeated from about 60N to the equator. It all comes from the same set
of data, processed at the same time with the script above and served
by one single layer in mapserver, so I can't make head or tails of it.

All ideas will be greatly appreciated.

Z


-- 
The best defence against logic is ignorance. The next best
is stupidity. Both can be used simultaneously.