[GRASS-user] r.stats - wrong area values from MOLA dataset

[Carlos Grohmann]

Hello all

I'm trying to get some stats for Mars using a MOLA dataset.

Since GRASS 7.0beta is having some problems with the most recent matplotlib
i added the planetary spheroids from v.7.0 into v.6.4.4 ellipse.table file.
Then I created a Location using Equirectangular (eqc) projection, no datum
definition (dx=dy=dz=0) and mars (IAU2000) spheroid.

After that I run, through python, r.stats and make a list of the area
values:

gdem = 'mola'
bwidth = 10 #10m bins for histogram

rinfo = grass.parse_command('r.info', map=gdem, flags='r')
g_min = float(rinfo['min'])
g_max = float(rinfo['max'])

grass.run_command('g.region', rast=gdem, flags='pa')
nbins = round((g_max - g_min)/bwidth) # bin size
mapstats = grass.read_command('r.stats', input = gdem, \
    fs = ',', nsteps = nbins, flags = 'an')

histlist = mapstats.splitlines()
elevlist = []
arealist = []
for pair in histlist:
    try:
        pairlist = pair.split(',')
        elev = int(float(pairlist[0]))
        area = float(pair.split(',')[1])
        elevlist.append(elev)
        arealist.append(area)
    except:
        pass

np.sum(arealist)

this last step will give me:
>>> np.sum(arealist)
64793.652361588101

Of course this is wrong. Mars has a surface area of about 144.8E6km2. I
know that my steps are working because I get correct values for Earth:

#with ETOPO5 (r.stats returns values in m2)
>>> np.sum(arealist)/10e6
507756893328616.19


So, what am I missing?

thanks

Carlos



-- 
Prof. Carlos Henrique Grohmann
Institute of Energy and Environment - Univ. of São Paulo, Brazil
- Digital Terrain Analysis | GIS | Remote Sensing -

http://carlosgrohmann.com
http://orcid.org/0000-0001-5073-5572
________________
Can't stop the signal.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20140701/8f49169c/attachment-0001.html>