[GRASS-user] Issues with resampling to equal pixel size for FCELL rasters

maning sambale emmanuel.sambale at gmail.com
Mon Dec 22 22:46:48 PST 2014


Hi,

I have a DTM I want to resample. The original data was in lat/lon SRS
which I re-projected to UTM and imported to GRASS GIS.

GRASS 6.4.3 (pampanga_wshed):~ > r.info ifsar_dtm
 +----------------------------------------------------------------------------+
 | Layer:    ifsar_dtm                      Date: Mon Dec 22 10:07:06 2014    |
 | Mapset:   PERMANENT                      Login of Creator: maning          |
 | Location: pampanga_wshed                                                   |
 | DataBase: /Users/maning/GRASSDATA/                                         |
 | Title:     ( ifsar_dtm )                                                   |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 255             |
 |   Data Type:    FCELL                                                      |
 |   Rows:         44440                                                      |
 |   Columns:      53950                                                      |
 |   Total Cells:  2397538000                                                 |
 |        Projection: UTM (zone 51)                                           |
 |            N: 1842577.192205    S: 1620332.25097015   Res: 5.00101128      |
 |            E: 420942.75364933    W: 151243.07137293   Res: 4.99906733      |
 |   Range of data:    min = -10000  max = 3.402823e+38                       |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.in.gdal                                                  |
 |                                                                            |
 |   Comments:                                                                |
 |    r.in.gdal -e input="ifsar_dtm_utmz51n.vrt" output="ifsar_dtm"           |
 |                                                                            |
 +----------------------------------------------------------------------------+

Because of the origin SRS, the pixel size is not the same for ew_res and nw_res.
I tried running several resampling techniques to match the pixel size
equal to 5.

First, I tried r.resamp.interp.

g.region rast=ifsar_dtm_subset -pm
projection: 1 (UTM)
zone:       51
datum:      wgs84
ellipsoid:  wgs84
north:      1683792
south:      1645429
west:       220827
east:       268500
nsres:      5.00104289
ewres:      4.99926594
rows:       7671
cols:       9536
cells:      73150656

g.region res=5 -ap
projection: 1 (UTM)
zone:       51
datum:      wgs84
ellipsoid:  wgs84
north:      1683795
south:      1645425
west:       220825
east:       268500
nsres:      5
ewres:      5
rows:       7674
cols:       9535
cells:      73171590

r.resamp.interp input=ifsar_dtm_subset
output=ifsar_dtm_subset_interp_bic method=bicubic --overwrite
r.info ifsar_dtm_subset_interp_bic
 +----------------------------------------------------------------------------+
 | Layer:    ifsar_dtm_subset_interp_bic    Date: Tue Dec 23 13:17:51 2014    |
 | Mapset:   test1                          Login of Creator: maning          |
 | Location: pampanga_wshed                                                   |
 | DataBase: /Users/maning/GRASSDATA/                                         |
 | Title:    Resample by bicubic interpolation ( ifsar_dtm_subset_interp_bic  |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 255             |
 |   Data Type:    DCELL                                                      |
 |   Rows:         7674                                                       |
 |   Columns:      9535                                                       |
 |   Total Cells:  73171590                                                   |
 |        Projection: UTM (zone 51)                                           |
 |            N:    1683795    S:    1645425   Res:     5                     |
 |            E:     268500    W:     220825   Res:     5                     |
 |   Range of data:    min = -2.72093853362067  max = 1157.05376343899        |
 |                                                                            |
 |   Data Source:                                                             |
 |    ifsar_dtm_subset                                                        |
 |    Source map NS res: 5.00104289   EW res: 4.99926594                      |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.resamp.interp                                            |
 |                                                                            |
 |   Comments:                                                                |
 |    r.resamp.interp input="ifsar_dtm_subset" output="ifsar_dtm_subset_in\   |
 |    terp_bic" method="bicubic"                                              |
 |                                                                            |
 +----------------------------------------------------------------------------+


The output resolution matches.

Then, I tried r.resamp.rst.

g.region rast=ifsar_dtm_subset -pm
projection: 1 (UTM)
zone:       51
datum:      wgs84
ellipsoid:  wgs84
north:      1683792
south:      1645429
west:       220827
east:       268500
nsres:      5.00104289
ewres:      4.99926594
rows:       7671
cols:       9536
cells:      73150656

r.resamp.rst input=ifsar_dtm_subset elev=ifsar_dtm_subset_rrst
ew_res=5 ns_res=5 --overwrite
r.info ifsar_dtm_subset_rrst
 +----------------------------------------------------------------------------+
 | Layer:    ifsar_dtm_subset_rrst          Date: Tue Dec 23 14:13:49 2014    |
 | Mapset:   test1                          Login of Creator: maning          |
 | Location: pampanga_wshed                                                   |
 | DataBase: /Users/maning/GRASSDATA/                                         |
 | Title:     ( ifsar_dtm_subset_rrst )                                       |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 255             |
 |   Data Type:    FCELL                                                      |
 |   Rows:         7673                                                       |
 |   Columns:      9535                                                       |
 |   Total Cells:  73162055                                                   |
 |        Projection: UTM (zone 51)                                           |
 |            N:    1683792    S:    1645429   Res: 4.99973935                |
 |            E:     268500    W:     220827   Res: 4.99979025                |
 |   Range of data:    min = -2.759012  max = 1157.082                        |
 |                                                                            |
 |   Data Source:                                                             |
 |    raster map ifsar_dtm_subset                                             |
 |                                                                            |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.resamp.rst                                               |
 |                                                                            |
 |   Comments:                                                                |
 |    tension=617.194349                                                      |
 |    dnorm=64.809407, zmult=1.000000                                         |
 |    KMAX=50, KMIN=35, errtotal=0.003577                                     |
 |    zmin_data=-2.600000, zmax_data=1157.017822                              |
 |    zmin_int=-2.759012, zmax_int=1157.081551                                |
 |                                                                            |
 +----------------------------------------------------------------------------+

Really close but still doesn't match.  To visualize the output I
created aspect layers for each respection resampling method, but both
layer still show checkerboard effects.

https://drive.google.com/file/d/0B-2WZQ1DwK_xbUpBTFJzcFdXV2M/view?usp=sharing
Any ideas on what I missed?

BTW, off topic, is there a way to limit the precision to say 3 decimal
places?  Centimeter vertical unit should be more than enough for my
intended use.


-- 
cheers,
maning
------------------------------------------------------
"Freedom is still the most radical idea of all" -N.Branden
blog: http://epsg4253.wordpress.com/
------------------------------------------------------


More information about the grass-user mailing list