[GRASS-user] Issues with resampling to equal pixel size for FCELL rasters
Markus Metz
markus.metz.giswork at gmail.com
Tue Dec 23 06:07:33 PST 2014
On Tue, Dec 23, 2014 at 7:46 AM, maning sambale
<emmanuel.sambale at 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 at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user
More information about the grass-user
mailing list