Issues when reprojecting from a GCS to a PCS

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Issues when reprojecting from a GCS to a PCS

Hanlie Pretorius
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
_______________________________________________
grass-user mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-user
Reply | Threaded
Open this post in threaded view
|

Re: Issues when reprojecting from a GCS to a PCS

Neil Best
Hanlie, judging from a quick scan of your problem you might consider using gdalbuildvrt to build a virtual mosaic of the inputs and gdalwarp to reproject the result using the VRT driver.  I have found this to be quite elegant and efficient since you have accomplished your goal by creating two tiny XML files, although you may have to edit them by hand if you want to move the VRT files and adjust the paths to the real data, but it can be done.  You can even open the VRTs in QGIS if you want to visualize before importing.  In my experience the GDAL utilities set the resolution of in a more sensible way that is described clearly in the docs.  If this is a one-time operation then there is no down side to running it through two layers of VRT vs. writing the reprojected data out to disk just to have r.in.gdal write it yet again to your GRASS database.  Let me know if that's enough to get you going or not.  Have fun!


Neil