[GRASS-dev] Reproject, change resolution and patch: doubts....

GRASS GIS GRASS GIS developer
1

Hi all,

I would like to submit to you the following “work plan” that I figured out to do a particular job:

I need to apply the module r.horizon (add-on from JRC) that computes the angular height of terrain horizon from a given point and a given DEM raster.
That DEM should be as large as possible to let the module create an horizon profile as much realistic as possible.
I have two types of GIS data of the interested area to work on:

A) the SRTM database, 90m resolution, WGS84
B) a local GIS database, with a resolution of 20m, Gauss-Boaga Rome40 West (preferred, because of the better resoultion)

all the data are in arc-info ASCII format.
Since I don’t want to create confusion within different data and projects, I decided to create two main locations: a Gass-Boaga local for the local database, in which import ASCII data usign r.in.arc module, and an WGS84 location for the SRTM database (using r.in.arc here too).

While the SRTM data are organanized in big maps, where a single map covers all the interested area, the Local database is fragmented in multiple small (but contiguos) maps; because of that I need to patch all of them to have a sufficient working area. But that’s not still enough, since some points of interest lay on the borders of such local map: that’s why I need to patch the local database big map (all the maps patched together) with the SRTM map (only one map is needed), with the following goals:

  1. obtain a map with the local map data, where those are present, and with the SRTM data where the local data are not present (null)
  2. do not alter the original data of both locations during patching or resampling (for example for SRTM passing from 90m resolution to 20m)

I try to briefly summarize the steps of the procedure I figured out to reach that goals:

  1. Create a Gauss-Boaga location (call it GB location) and import the local database (native resolution 20m) with r.in.arc
  2. Set the region in order to contain all the GB maps and patch all the imported maps (contiguous) to create the big map
  3. Create a WGS84 location (call it SRTM location) and import the needed DRTM map (native resolution 90m) with r.in.arc
  4. Create another WGS84 location (all it JOB location; I create it new because I prefer to preserve the original locations) and set the region to fit the interested area
  5. In the JOB Location: r.proj the GB big map, maintaining the orginal resolution (20m)
  6. r.proj the SRTM map, setting the resolution to 20m (against the original of 90m) during the reprojection (to let me patch it with the GB map that has a 20m res)
  7. r.patch the GB big map and the SRTM map, putting the GB raster as first in the input parameter to achieve the goal n°1

My doubt is: moving from the 90m resolution to 20m in the SRTM map I will oversample, right? but I will actually maintain the original data (that it is what I want).
Instead, if I resample the SRTM data (using one of r.resamp.* modules) I would obtain a smooth DEM, but I would also alterate the original data, introducing dangerous fake data, and so possible critical errors.

Do you agree?

Thanks for all,

Marco

2

Hi Helena,

Thanks for your reply.

just a few notes. If you set the resolution to 20m, GRASS will resample

your 90m SRTM to 20m using nearest neighbor, creating a DEM with "steps" -
something line this (a) is automatic nn resampling, b) is reinterpolation by
rst) http://www.grassbook.org/gallery/chapters/s050303f010_resamp.jpg

Yes, yesterday I tried the described procedure and I noticed the same result
NVIZing the reprojected dem. BTW I don't need to produce 3D maps, but I need
to introduce the less data distorsion as possible. The NN method should do
the job for me. I'm not sure that an interpolation would be a good solution,
but I talk with a strict mathematic/electronic approach, that's my
professional training; probably, talking about topography, interpolations
rarely produce remarkable errors.

This said, I prefer the 20m data, against the SRTM ones, because I analyzed
a very well known (small) area with NVIZ, using both SRTM and 20m data, and
I noticed relevant errors with SRTM data (I cannot pretend much from it,
it's a SAR shuttle mission... it's actually a miracle that we have it, that
it works pretty fine and that is free).

Summarizing:
1) r.proj method=nearest
2) r.resamp.rst

Right?

Thanks,

Marco

-----Messaggio originale-----
Da: Helena Mitasova [mailto:hmitaso@unity.ncsu.edu]
Inviato: mercoledì 7 maggio 2008 18.25
A: marco.pasetti@alice.it
Oggetto: Re: [GRASS-dev] Reproject, change resolution and patch: doubts....

Marko,

just a few notes. If you set the resolution to 20m, GRASS will resample your
90m SRTM to 20m using nearest neighbor, creating a DEM with "steps" -
something line this (a) is automatic nn resampling, b) is reinterpolation by
rst) http://www.grassbook.org/gallery/chapters/s050303f010_resamp.jpg
Dylan has a nice comparison of effects of different resampling and
reinterpolation methods on his web site, if you want to learn more.

If you have a GRASS book and the nc_spm data set - it has both SRTM data and
lidar based 10m and 1m resolution DEMs so you can see how SRTM works, there
are several examples in the book - e.g. SRTM is consistently higher and I
did viewshed from 30m SRTM and lidar-based 30m DEM - same resolution but the
results were quite different. You need to keep in mind that SRTM maps the
surface with vegetation on it while your 20m DEM is probably bare earth
surface. Also look at the metadata for both to compare their vertical
accuracy, some smoothing of SRTM may actually be justified (unless you want
the vegetation and have some means to add it to your 20m DEM). When you
decide to reinterpolate you can compare your re-interpolated DEM with the
original data (e.g. the 90m grid centers) and make sure that any smoothing
does not exceed the accuracy of your original data.

Just for illustration: here is a stream network derived from 10m resolution
IFSARE patched with 90m SRTM where IFSARE data were missing, all
reinterpolated to 30m resolution to get the river flow seamlessly
back-and-forth along the patch line.
http://skagit.meas.ncsu.edu/~helena/measwork/panama/pacora30maccum.jpg
some SRTM data can be pretty noisy:
http://skagit.meas.ncsu.edu/~helena/measwork/nrc/mumbai_srtm.png

Helena

Do you agree?

Thanks for all,

Marco
_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

3

Hi Helena,

  1. r.proj method=nearest

good, the basic step has been fixed :slight_smile:

  1. r.resamp.rst
  • you don’t need this step if the default nearest neighbor works for
    your application
  • it is useful only if you needed smooth connection between the 20m
    DEM and SRTM. ( r.resamp.rst does interpolation using a spline
    function - you can make it pass exactly through the given data points
    by setting smoothing =0.)

Actually I already thought that, in any case, I should at least fix the interferences among the two different maps. That means that I should at least apply r.resamp.rst along the line where the two maps collide. Having the 20m raster called, for example, 20mDEM, I figured out (and already tested) how to create a polygon matching the exact boudaries of that map (exactly an area vector layer matching the whole area covered by the 20mDEM map):

create a raster from 20mDEM with NULL if NULL and 0 if >=0
transform that raster into a vector layer faeaturing the area covered by it

r.mapcalc 20mDEM_surface=‘if(20mDEM,0)’
r.to.vect feauture=area

Then I could apply r.resamp.rst on a buffer along the polyline that describe the border of the above area… but I’m not sure how to do that, I still need to carefully read the online manuals (unfortunately I still don’t have the book of Markus; GFOSS.it association kindly gave me a copy… but I still have to get it…)

Thanks,

Marco

Da: Helena Mitasova [mailto:hmitaso@unity.ncsu.edu]
Inviato: gio 08/05/2008 2.46
A: marco.pasetti@alice.it
Oggetto: Re: R: [GRASS-dev] Reproject, change resolution and patch: doubts…

On May 7, 2008, at 3:35 PM, Marco Pasetti wrote:

Hi Helena,

Thanks for your reply.

just a few notes. If you set the resolution to 20m, GRASS will
resample
your 90m SRTM to 20m using nearest neighbor, creating a DEM with
“steps” -
something line this (a) is automatic nn resampling, b) is
reinterpolation by
rst) http://www.grassbook.org/gallery/chapters/s050303f010_resamp.jpg

Yes, yesterday I tried the described procedure and I noticed the
same result
NVIZing the reprojected dem. BTW I don’t need to produce 3D maps,
but I need
to introduce the less data distorsion as possible. The NN method
should do
the job for me. I’m not sure that an interpolation would be a good
solution,
but I talk with a strict mathematic/electronic approach, that’s my
professional training; probably, talking about topography,
interpolations
rarely produce remarkable errors.

This said, I prefer the 20m data, against the SRTM ones, because I
analyzed
a very well known (small) area with NVIZ, using both SRTM and 20m
data, and
I noticed relevant errors with SRTM data (I cannot pretend much
from it,
it’s a SAR shuttle mission… it’s actually a miracle that we have
it, that
it works pretty fine and that is free).

Summarizing:

  1. r.proj method=nearest
  1. r.resamp.rst
  • you don’t need this step if the default nearest neighbor works for
    your application
  • it is useful only if you needed smooth connection between the 20m
    DEM and SRTM. ( r.resamp.rst does interpolation using a spline
    function - you can make it pass exactly through the given data points
    by setting smoothing =0.)

Helena

Right?

Thanks,

Marco

-----Messaggio originale-----
Da: Helena Mitasova [mailto:hmitaso@unity.ncsu.edu]
Inviato: mercoledì 7 maggio 2008 18.25
A: marco.pasetti@alice.it
Oggetto: Re: [GRASS-dev] Reproject, change resolution and patch:
doubts…

Marko,

just a few notes. If you set the resolution to 20m, GRASS will
resample your
90m SRTM to 20m using nearest neighbor, creating a DEM with “steps” -
something line this (a) is automatic nn resampling, b) is
reinterpolation by
rst) http://www.grassbook.org/gallery/chapters/s050303f010_resamp.jpg
Dylan has a nice comparison of effects of different resampling and
reinterpolation methods on his web site, if you want to learn more.

If you have a GRASS book and the nc_spm data set - it has both SRTM
data and
lidar based 10m and 1m resolution DEMs so you can see how SRTM
works, there
are several examples in the book - e.g. SRTM is consistently higher
and I
did viewshed from 30m SRTM and lidar-based 30m DEM - same
resolution but the
results were quite different. You need to keep in mind that SRTM
maps the
surface with vegetation on it while your 20m DEM is probably bare
earth
surface. Also look at the metadata for both to compare their vertical
accuracy, some smoothing of SRTM may actually be justified (unless
you want
the vegetation and have some means to add it to your 20m DEM). When
you
decide to reinterpolate you can compare your re-interpolated DEM
with the
original data (e.g. the 90m grid centers) and make sure that any
smoothing
does not exceed the accuracy of your original data.

Just for illustration: here is a stream network derived from 10m
resolution
IFSARE patched with 90m SRTM where IFSARE data were missing, all
reinterpolated to 30m resolution to get the river flow seamlessly
back-and-forth along the patch line.
http://skagit.meas.ncsu.edu/~helena/measwork/panama/pacora30maccum.jpg
some SRTM data can be pretty noisy:
http://skagit.meas.ncsu.edu/~helena/measwork/nrc/mumbai_srtm.png

Helena

Do you agree?

Thanks for all,

Marco

grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

4

Marco Pasetti pisze:
Summarizing:

1) r.proj method=nearest

> 2) r.resamp.rst

Right?

Marco,

If you care about the accuracy very much:

# in ll/WGS84 location:
g.region rast=your_srtm_tile
r.to.vect -z in=your_srtm_tile feature=point out=your_srtm_tile_vect

# in your projected location
v.proj
g.region vect=your_srtm_tile_vect res=20 -a
v.surf.rst (or else)

Raster to vector followed by reprojecting the vector preserves the
original SRTM values perfectly in regard to their values and location
(given that you set the initial region correct). Also it lets you avoid
the distortion which takes place when reprojecting raster, due to
different grid spacing and skew of the source latlong and target
projected location.

Maciek

P.S.

Just for the record - it looks more like an issue to dicsuss on the
user ML.

5

Marco:

Actually I already thought that, in any case, I should at least fix the
interferences among the two different maps. That means that I should at
least apply r.resamp.rst along the line where the two maps collide.

Using r.series to do the patching might help smooth over any overlaps.

Hamish

      ____________________________________________________________________________________
Be a better friend, newshound, and
know-it-all with Yahoo! Mobile. Try it now. http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ