Hi Everyone,
I'm working in GRASS 6.4RC6 on Win XP and am trying to reproject TRMM
rainfall data from its native GCS (WGS84) to a local PCS. I'm using
the PCS because my other important data, specifically DEMs, are
provided in it. TRMM data come in 0.25°x0.25° squares.
My workfow is as follows:
1. Import the TRMM data into a WGS84 location using r.in.xyz. This
gives a raster with the following information:
+----------------------------------------------------------------------------+
| Layer: 3B42.000215.0@PERMANENT Date: Thu Jul 15 14:42:27 2010
| Mapset: PERMANENT Login of Creator: hanlie
| Location: world_wgs84
| DataBase: F:\Hanlie\grassdata
| Title: Raw x,y,z data binned into a raster grid by cell mean ( 3B42.000
| Timestamp: none
|---------------------------------------------------------------------------
|
| Type of Map: raster Number of Categories: 255
| Data Type: FCELL
| Rows: 6
| Columns: 3
| Total Cells: 18
| Projection: Latitude-Longitude
| N: 27:15S S: 28:45S Res: 0:15
| E: 28:45E W: 28E Res: 0:15
| Range of data: min = 0.000000 max = 0.070614
|
| Data Source:
| F:\Hanlie\UCT\M.Sc\Data\TRMM\2000\02_Februarie\3B42.000215.0.6.nc.lieb.
|
|
| Data Description:
| generated by r.in.xyz
|
| Comments:
| r.in.xyz input="F:\Hanlie\UCT\M.Sc\Data\TRMM\2000\02_Februarie\3B42.\
| 000215.0.6.nc.lieb.txt" output="3B42.000215.0" method="mean" type="F\
| CELL" x=2 y=1 z=3 zscale=1.0 percent=100
|
+----------------------------------------------------------------------------+
2. Set the region to this raster file and create a vector map from it
using v.in.region. This gives a vector with the following information:
+----------------------------------------------------------------------------+
| Layer: lieb_trmm_region@PERMANENT
| Mapset: PERMANENT
| Location: world_wgs84
| Database: F:\Hanlie\grassdata
| Title:
| Map scale: 1:1
| Map format: native
| Name of creator: hanlie
| Organization:
| Source date: Thu Jul 15 14:55:13 2010
|----------------------------------------------------------------------------|
| Type of Map: vector (level: 2)
|
| Number of points: 0 Number of areas: 1
| Number of lines: 0 Number of islands: 1
| Number of boundaries: 1 Number of faces: 0
| Number of centroids: 1 Number of kernels: 0
|
| Map is 3D: No
| Number of dblinks: 0
|
| Projection: Lat/Lon
| N: 27:15S S: 28:45S
| E: 28:45E W: 28E
|
| Digitization threshold: 0
| Comments:
|
+----------------------------------------------------------------------------+
3. Change to the PCS location and reproject the TRMM region using
v.proj. This gives a vector with the following information:
+----------------------------------------------------------------------------+
| Layer: lieb_trmm_region@PERMANENT
| Mapset: PERMANENT
| Location: sa_tm_29deg_E
| Database: F:\Hanlie\grassdata
| Title:
| Map scale: 1:1
| Map format: native
| Name of creator: hanlie
| Organization:
| Source date: Thu Jul 15 14:55:13 2010
|----------------------------------------------------------------------------|
| Type of Map: vector (level: 2)
|
| Number of points: 0 Number of areas: 1
| Number of lines: 0 Number of islands: 1
| Number of boundaries: 1 Number of faces: 0
| Number of centroids: 1 Number of kernels: 0
|
| Map is 3D: No
| Number of dblinks: 0
|
| Projection: Transverse Mercator
| N: -3015356.39959204 S: -3181970.90748457
| E: -24418.1519951 W: -99037.39660484
|
| Digitization threshold: 0
| Comments:
|
+----------------------------------------------------------------------------+
4. Set the region of the PCS to the TRMM region:
g.region -p vect=lieb_trmm_region@PERMANENT nsres=25 ewres=25
which gives:
projection: 99 (Transverse Mercator)
zone: 0
datum: ** unknown (default: WGS84) **
ellipsoid: wgs84
north: -3015356.39959204
south: -3181970.90748457
west: -99037.39660484
east: -24418.1519951
nsres: 24.99842579
ewres: 24.9980719
rows: 6665
cols: 2985
cells: 19895025
I use this resolution because I don't really know what else to use.
When I reproject the TRMM GRID into the PCS, I get polygons with
differing areas, so I can't use the TRMM grid size as resolution:
v.report map=lieb_trmm_squares@PERMANENT option=area units=kilometers
Displaying column types/names for database connection of layer 1:
cat|LabelX|LabelY|LblOffsetX|LblOffsetY|Label|gRow|gColumn|RowCol|area
1|28|-28.5|-20|-20|28||||679.062503849138
2|28|-28.25|-20|-20|28||||680.623828116383
3|28|-28|-20|-20|28||||682.172234597964
4|28|-27.75|-20|-20|28||||683.70769882544
5|28.25|-28.75|-20|-20|28.25||||677.428350142518
6|28.25|-28.5|-20|-20|28.25||||679.002141385312
7|28.25|-28.25|-20|-20|28.25||||680.563040996854
8|28.25|-28|-20|-20|28.25||||682.111024322052
9|28.25|-27.75|-20|-20|28.25||||683.646066944101
10|28.25|-27.5|-20|-20|28.25||||685.168144684034
11|28.5|-28.5|-20|-20|28.5||||678.961903132767
12|28.5|-28.25|-20|-20|28.5||||680.522519692921
13|28.5|-28|-20|-20|28.5||||682.070220966439
14|28.5|-27.75|-20|-20|28.5||||683.604982570918
15|28.5|-27.5|-20|-20|28.5||||685.126780361735
5. Reproject the TRMM rainfall from the GCS to the PCS, which gives no
data in the file. The extents also seem to be wrong:
+----------------------------------------------------------------------------+
| Layer: 3B42.000215.0@PERMANENT Date: Thu Jul 15 15:01:05 2010
| Mapset: PERMANENT Login of Creator: hanlie
| Location: sa_tm_29deg_E
| DataBase: F:\Hanlie\grassdata
| Title: ( 3B42.000215.0 )
| Timestamp: none
|----------------------------------------------------------------------------|
|
| Type of Map: raster Number of Categories: 255
| Data Type: FCELL
| Rows: 5552
| Columns: 1996
| Total Cells: 11081792
| Projection: Transverse Mercator
| N: -3029230.52590567 S: -3168021.78589357 Res: 24.99842579
| E: -36667.20722417 W: -86563.35872871 Res: 24.9980719
| Range of data: min = nan max = nan
|
| Data Description:
| generated by r.proj
|
| Comments:
| r.proj input="3B42.000215.0" location="world_wgs84" mapset="PERMANEN\
| T" method="cubic" resolution=25.0
|
+----------------------------------------------------------------------------+
Can someone perhaps see where I'm going wrong? I have hundreds of
these files to import, so I need to get the workflow right and then
automate the task.
Thanks
Hanlie