[GRASS-user] r.report units in LatLong location

[Hamish]

> Gordon Keith wrote:
> > I'm running r.report and specifying units=me in a Lat-Long location but the
> > values it returns appear to be in square degrees not square meters.
> > 
> > Is there a way to get it to report in square meters?

I have now tried it, it seems ok for sq miles. r.buffer seems LL safe too.


# etopo2 dataset:
g.region rast=etopo2

# create land/sea reclass map
r.reclass in=etopo2 out=landsea << EOF
-30000 thru -1 = 0  sea
0 thru 30000   = 1 land
EOF

# zoom to NZ
g.region n=-25 s=-60 w=155 e=190

# crop out NZ
r.mapcalc nzland='if(landsea, 1, null())'

# create a 100naut mile buffer  (200nm EEZ has messy overlaps)
r.buffer in=nzland out=nzland_buf100nm unit=nautmiles dist=100

# looks good, northern islands are buffered by circles,
#  southern islands buffered by ovals

# d.measure checks out ok



# zoom in on oval-looking Macquarie Isl. buffer
g.region n=-52 s=-57 w=155 e=163

r.report -nh nzland_buf100nm units=c,k,mi
+-------------------------------------------------------------------------+
|              Category Information           | cell|  square  |    square|
|#|description                                |count|kilometers|     miles|
|-------------------------------------------------------------------------|
|1|distances calculated from these locations. |    3|     24.00|     9.265|
|2|0-100 nautmiles. . . . . . . . . . . . . . |13977|111,778.53|43,157.689|
|-------------------------------------------------------------------------|
|TOTAL                                        |13980|111,802.52|43,166.954|
+-------------------------------------------------------------------------+


# zoom in on circular Lord Howe Isl. buffer
g.region n=-29 s=-34 w=156 e=162

r.report -nh nzland_buf100nm units=c,k,mi
+-------------------------------------------------------------------------+
|              Category Information           | cell|  square  |    square|
|#|description                                |count|kilometers|     miles|
|-------------------------------------------------------------------------|
|1|distances calculated from these locations. |    1|     11.70|     4.517|
|2|0-100 nautmiles. . . . . . . . . . . . . . | 9210|107,722.60|41,591.698|
|-------------------------------------------------------------------------|
|TOTAL                                        | 9211|107,734.30|41,596.214|
+-------------------------------------------------------------------------+

$ units '41596.214 miles^2' 'km^2'
        * 107733.7

hmph. that's not exactly the same as r.report, but close.

A small tweak in CVS to make things more precise, and all becomes even:

--- prt_report.c        30 Sep 2007 07:08:23 -0000      2.4
+++ prt_report.c        16 Nov 2007 04:52:14 -0000
@@ -63,7 +63,7 @@
        case ACRES:
            unit[i].label[0] = "";
            unit[i].label[1] = "acres";
-           unit[i].factor   = 2.471e-4;
+           unit[i].factor   = 2.4710439e-4;
            break;
 
        case HECTARES:
@@ -75,7 +75,7 @@
        case SQ_MILES:
            unit[i].label[0] = "square";
            unit[i].label[1] = " miles";
-           unit[i].factor   = 3.861e-7;
+           unit[i].factor   = 3.8610216e-7;
            break;
 
        default:



Just by hand for a 100nm circle,

100nm area in sq mi. is
  PI*(1.1507794mi/nm * 100)^2  =  41,604mi^2

The above examples are around a few 2nm cells, so the radius is a little
bigger:
  PI*(1.1507794 * 102)^2  =  43,284 mi^2


> If it never worked it will be a case of calculating the cell area per row for
> LL locations [nsres*60nm->m * (ewres*60nm->m * cos(lat))]. I don't really
> feel qualified to anticipate what sort of systematic cumulative errors that
> might introduce or what a better equal-area projection to use might be on
> such a continental scale.

I though about a new libgis function-
  double area = G_area_of_ll_cell(lat, ewres, nsres);
which would calculate the area of a trapezoid approximating a single LL cell in
sq. meters. That could be calculated per LL raster row. That would be better if
cell size was >1 deg or so when you're assuming the width at the center of the
cell is the average of the width at the top and bottom of the cell. Well, it's
not linear so that's really not much of an improvement on the current. I guess
you could split up the cell and calc area for the top half and bottom half from
two trapezoids.. of course you can take this further and make more and more
polygons but I guess it would be best to forget that way and just do the area
calc properly using spherical math. shrug




Anyway, AFAICT the r.report code looks ok and it's beer o'clock over here..


Hamish



      ____________________________________________________________________________________
Be a better sports nut!  Let your teams follow you 
with Yahoo Mobile. Try it now.  http://mobile.yahoo.com/sports;_ylt=At9_qDKvtAbMuh1G1SQtBI7ntAcJ