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

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/
------------------------------------------------------

On Tue, Dec 23, 2014 at 7:46 AM, maning sambale
<emmanuel.sambale@gmail.com> wrote:

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.

You have probably used some automated or on-the-fly reprojection
method. As you noticed, the results are not satisfying.

The simplest solution would be to use gdalwarp or GRASS r.proj to do
the reprojection, specifying extents and resolution for gdalwarp or
setting the region first for r.proj. Also take care of the resampling
method with gdalwarp and r.proj. Default is nearest neighbor which you
probably don't want for a DEM.

Markus M

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/
------------------------------------------------------
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user