[GRASS-user] r.proj issue

GRASS GRASS user
1

Reprojecting a 10m DEM fails. Following are the command and input/output
projection and map description information. Screenshots of the input and
output maps are attached.

I need to understand why the output is wrong.

Command:
r.proj loc=Oregon map=topography in=n46w123 met=lanczos_f mem=16000 --o

Input location:
name: unnamed
ellps: grs80
proj: lcc
lat_1: 46
lat_2: 44.33333333333334
lat_0: 43.66666666666666
lon_0: -120.5
x_0: 2500000
y_0: 0
no_defs: defined
unit: Meter
units: Meters
meters: 1

Cols: 502 (80961)
Rows: 39234 (113148)
North: 262349.035904 (262349.035904)
South: 223115.062541 (149201.112724)
West: 2302891.866970 (2302891.866970)
East: 2303393.867032 (2383852.876897)
EW-res: 1.000000
NS-res: 0.999999

Output location:
name: NAD83(HARN) / Oregon North
datum: nad83harn
ellps: grs80
proj: lcc
lat_1: 46
lat_2: 44.33333333333334
lat_0: 43.66666666666666
lon_0: -120.5
x_0: 2500000
y_0: 0
no_defs: defined
unit: meter
units: meters
meters: 1

Cols: 500 (71702)
Rows: 39232 (50757)
North: 262348.951018 (273873.984569)
South: 223116.836808 (223116.836808)
West: 2302892.017711 (2231689.553588)
East: 2303392.020970 (2303392.020970)
EW-res: 1.000007
NS-res: 1.000003Cols: 500 (71702)
Rows: 39232 (50757)
North: 262348.951018 (273873.984569)
South: 223116.836808 (223116.836808)
West: 2302892.017711 (2231689.553588)
East: 2303392.020970 (2303392.020970)
EW-res: 1.000007
NS-res: 1.000003

Rich

(attachments)

input-dem.png
output-dem.png

2

Rich Shepard wrote

Reprojecting a 10m DEM fails. Following are the command and input/output
projection and map description information. Screenshots of the input and
output maps are attached.

I need to understand why the output is wrong.

Command:
r.proj loc=Oregon map=topography in=n46w123 met=lanczos_f mem=16000 --o

Input location:
name: unnamed
ellps: grs80
proj: lcc
lat_1: 46
lat_2: 44.33333333333334
lat_0: 43.66666666666666
lon_0: -120.5
x_0: 2500000
y_0: 0
no_defs: defined
unit: Meter
units: Meters
meters: 1

Cols: 502 (80961)
Rows: 39234 (113148)
North: 262349.035904 (262349.035904)
South: 223115.062541 (149201.112724)
West: 2302891.866970 (2302891.866970)
East: 2303393.867032 (2383852.876897)
EW-res: 1.000000
NS-res: 0.999999

Output location:
name: NAD83(HARN) / Oregon North
datum: nad83harn
ellps: grs80
proj: lcc
lat_1: 46
lat_2: 44.33333333333334
lat_0: 43.66666666666666
lon_0: -120.5
x_0: 2500000
y_0: 0
no_defs: defined
unit: meter
units: meters
meters: 1

Cols: 500 (71702)
Rows: 39232 (50757)
North: 262348.951018 (273873.984569)
South: 223116.836808 (223116.836808)
West: 2302892.017711 (2231689.553588)
East: 2303392.020970 (2303392.020970)
EW-res: 1.000007
NS-res: 1.000003Cols: 500 (71702)
Rows: 39232 (50757)
North: 262348.951018 (273873.984569)
South: 223116.836808 (223116.836808)
West: 2302892.017711 (2231689.553588)
East: 2303392.020970 (2303392.020970)
EW-res: 1.000007
NS-res: 1.000003

Rich
_______________________________________________
grass-user mailing list

grass-user@.osgeo

https://lists.osgeo.org/mailman/listinfo/grass-user

input-dem.png (296K)
<http://osgeo-org.1560.x6.nabble.com/attachment/5416106/0/input-dem.png>
output-dem.png (28K)
<http://osgeo-org.1560.x6.nabble.com/attachment/5416106/1/output-dem.png>

a good thing to see if reprojection of raster data works correctly, is
comparing before/after r.proj with gdalwarp (1) results.

(1) https://gdal.org/programs/gdalwarp.html#gdalwarp

-----
best regards
Helmut
--
Sent from: http://osgeo-org.1560.x6.nabble.com/Grass-Users-f3884509.html

3

On Fri, 13 Sep 2019, Helmut Kudrnovsky wrote:

a good thing to see if reprojection of raster data works correctly, is
comparing before/after r.proj with gdalwarp (1) results.
(1) https://gdal.org/programs/gdalwarp.html#gdalwarp

Helmut,

I'll try gdalwarp. Scanning that web page suggests it does what r.proj does.
What puzzled me is why one raster map correctly projected while the second
didn't. Perhaps I need to download the source for the second one again. I'll
report results.

These are large files so downloading and processing are rather slow with my
low-speed FiOS connection and fast desktop.

Thanks again,

Rich

4

On Sat, Sep 14, 2019 at 1:15 AM Helmut Kudrnovsky <hellik@web.de> wrote:

Rich Shepard wrote

Reprojecting a 10m DEM fails. Following are the command and input/output
projection and map description information. Screenshots of the input and
output maps are attached.

I need to understand why the output is wrong.

Command:
r.proj loc=Oregon map=topography in=n46w123 met=lanczos_f mem=16000 --o

Input location:
name: unnamed
ellps: grs80
proj: lcc
lat_1: 46
lat_2: 44.33333333333334
lat_0: 43.66666666666666
lon_0: -120.5
x_0: 2500000
y_0: 0
no_defs: defined
unit: Meter
units: Meters
meters: 1

Cols: 502 (80961)
Rows: 39234 (113148)
North: 262349.035904 (262349.035904)
South: 223115.062541 (149201.112724)
West: 2302891.866970 (2302891.866970)
East: 2303393.867032 (2383852.876897)
EW-res: 1.000000
NS-res: 0.999999

Output location:
name: NAD83(HARN) / Oregon North
datum: nad83harn
ellps: grs80
proj: lcc
lat_1: 46
lat_2: 44.33333333333334
lat_0: 43.66666666666666
lon_0: -120.5
x_0: 2500000
y_0: 0
no_defs: defined
unit: meter
units: meters
meters: 1

Cols: 500 (71702)
Rows: 39232 (50757)
North: 262348.951018 (273873.984569)
South: 223116.836808 (223116.836808)
West: 2302892.017711 (2231689.553588)
East: 2303392.020970 (2303392.020970)
EW-res: 1.000007
NS-res: 1.000003Cols: 500 (71702)
Rows: 39232 (50757)
North: 262348.951018 (273873.984569)
South: 223116.836808 (223116.836808)
West: 2302892.017711 (2231689.553588)
East: 2303392.020970 (2303392.020970)
EW-res: 1.000007
NS-res: 1.000003

Rich

grass-user mailing list

grass-user@.osgeo

https://lists.osgeo.org/mailman/listinfo/grass-user

input-dem.png (296K)
<http://osgeo-org.1560.x6.nabble.com/attachment/5416106/0/input-dem.png&gt;
output-dem.png (28K)
<http://osgeo-org.1560.x6.nabble.com/attachment/5416106/1/output-dem.png&gt;

a good thing to see if reprojection of raster data works correctly, is
comparing before/after r.proj with gdalwarp (1) results.

gdalwarp estimates the target extents unless you use the -te option. In contrast, r.proj uses the current region. r.proj can estimate the target extents with the -p or -g flag. I suggest to first check you current region in the target location/mapset and compare with r.proj -p or -g

Markus M

(1) https://gdal.org/programs/gdalwarp.html#gdalwarp

best regards
Helmut

Sent from: http://osgeo-org.1560.x6.nabble.com/Grass-Users-f3884509.html

grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

5

On Fri, 13 Sep 2019, Rich Shepard wrote:

Reprojecting a 10m DEM fails. Following are the command and input/output
projection and map description information. Screenshots of the input and
output maps are attached.

I need to understand why the output is wrong.

Stepping away from the issue for a while often lets in the light of
understanding, as it did in this case.

I ran r.info on the map in the project directory (the one with the bagette
for a region) and saw that the minimum and maximum elevation values were
zero. No wonder the region was so narrow and displayed nothing.

Stepping back one level to the reprojection of the imported raw data, r.info
showed that that that copy also had no elevation values.

The imported lon-lat data has a maximum elevation of 1614.5 feet (the Coast
Range of Oregon) and a minimum of 0 feet where the Pacific Ocean washes
ashore at high tide.

Re-projecting this lon-lat map to the Oregon state projection location
worked this second time. r.info confirms the maximum elevation value. Now
the system is crunching to reproject that map to the project location, and I
expect this one will also be good.

From now on I'll check r.info and v.info to troubleshoot unexpected map

issues before asking on this mail list.

Thanks, Helmut.

Regards,

Rich

6

On Sat, 14 Sep 2019, Rich Shepard wrote:

Stepping away from the issue for a while often lets in the light of
understanding, as it did in this case.

Wrong. The source and destination locations have the same projections:

Input:
name: NAD83(HARN) / Oregon North
datum: nad83harn
ellps: grs80
proj: lcc
lat_1: 46
lat_2: 44.33333333333334
lat_0: 43.66666666666666
lon_0: -120.5
x_0: 2500000
y_0: 0
no_defs: defined
towgs84: 0.000,0.000,0.000
unit: meter
units: meters
meters: 1

Output:
name: NAD83(HARN) / Oregon North
datum: nad83harn
ellps: grs80
proj: lcc
lat_1: 46
lat_2: 44.33333333333334
lat_0: 43.66666666666666
lon_0: -120.5
x_0: 2500000
y_0: 0
no_defs: defined
towgs84: 0.000,0.000,0.000
unit: meter
units: meters
meters: 1

The input map (the source for all statewide data) has elevations while the
output map (project-specific maps) does not (max=NULL). This suggests to me
that r.proj is the wrong tool because there's no reprojection required, only
copying from one location to another. But, the g.copy manual page says it
copies from one mapset to another mapset in the same location.

How do I copy maps from one location to another when the projection of both
is the same?

Regards,

Rich

7

Rich Shepard <rshepard@appl-ecosys.com> schrieb am So., 15. Sep. 2019, 01:47:

How do I copy maps from one location to another when the projection of both
is the same?

With r.pack/v.pack, then r.unpack/v.unpack.

Markus

8

On Sun, 15 Sep 2019, Markus Neteler wrote:

How do I copy maps from one location to another when the projection of both
is the same?

With r.pack/v.pack, then r.unpack/v.unpack.

Markus,

Well! I definitely missed recognizing these as the answer when I looked at
the list of modules. Perhaps because I've not before used the pack/unpack
pair. Thanks for the pointer.

Best regards,

Rich

9

On Sat, 14 Sep 2019, Markus Metz wrote:

gdalwarp estimates the target extents unless you use the -te option. In
contrast, r.proj uses the current region. r.proj can estimate the target
extents with the -p or -g flag. I suggest to first check you current
region in the target location/mapset and compare with r.proj -p or -g

Markus M, et al.:

I've worked on this for 2+ days and cannot get the DEM from the source
location to the target location when both have the same projection and
units.

Following the above still fails and I've no idea what to do now.

Source location:

r.proj -p loc=Oregon map=topography in=n45w121

Input map <n45w121@topography> in location <Oregon>:
Source cols: 80297
Source rows: 111356
Local north: 148340.8192959
Local south: 36984.78211177
Local west: 2459852.09770149
Local east: 2540149.13365527

Current (target) region:

g.region -p

projection: 99 (NAD83(HARN) / Oregon North)
zone: 0
datum: nad83harn
ellipsoid: grs80
north: 307220.87234928
south: -199204.29675669
west: 2161358.9640483
east: 2822158.25992759
nsres: 1.00000033
ewres: 1.00000045
rows: 506425
cols: 660799
cells: 334645133575

Reset current region to the source:

g.region n=148340.8192959 s=36984.78211177 w=2459852.09770149 e=2540149.13365527 cols=80297 rows=111356

Unpack the DEM in the current location/mapset;

r.unpack n45w121.pack

WARNING: Difference between PROJ_INFO file of packed map and of current
          location:
          name: NAD83(HARN) / Oregon North
          datum: nad83harn
          ellps: grs80
          proj: lcc
          lat_1: 46
          lat_2: 44.33333333333334
          lat_0: 43.66666666666666
          lon_0: -120.5
          x_0: 2500000
          y_0: 0
          no_defs: defined
          - towgs84: 0.000,0.000,0.000
ERROR: Projection of dataset does not appear to match current location. In
        case of no significant differences in the projection definitions,
        use the -o flag to ignore them and use current location definition.

If I correctly understand gdalwarp it functions like r.proj in reprojecting
a raster map from one location to another location. Here, both the source
and target locations have the same projection. Will gdalwarp work in this
case? Isn't r.pack/r.unpack the correct module pair?

Please advise me of what I should do now or what other information I can
provide.

Regards,

Rich

10

On Sun, Sep 15, 2019 at 8:39 PM Rich Shepard <rshepard@appl-ecosys.com> wrote:

On Sat, 14 Sep 2019, Markus Metz wrote:

gdalwarp estimates the target extents unless you use the -te option. In
contrast, r.proj uses the current region. r.proj can estimate the target
extents with the -p or -g flag. I suggest to first check you current
region in the target location/mapset and compare with r.proj -p or -g

Markus M, et al.:

I’ve worked on this for 2+ days and cannot get the DEM from the source
location to the target location when both have the same projection and
units.

Following the above still fails and I’ve no idea what to do now.

Source location:

r.proj -p loc=Oregon map=topography in=n45w121
Input map n45w121@topography in location :
Source cols: 80297
Source rows: 111356
Local north: 148340.8192959
Local south: 36984.78211177
Local west: 2459852.09770149
Local east: 2540149.13365527

Current (target) region:

g.region -p
projection: 99 (NAD83(HARN) / Oregon North)
zone: 0
datum: nad83harn
ellipsoid: grs80
north: 307220.87234928
south: -199204.29675669
west: 2161358.9640483
east: 2822158.25992759
nsres: 1.00000033
ewres: 1.00000045
rows: 506425
cols: 660799
cells: 334645133575

Reset current region to the source:

g.region n=148340.8192959 s=36984.78211177 w=2459852.09770149 e=2540149.13365527 cols=80297 rows=111356

You should now use
r.proj loc=Oregon map=topography in=n45w121

and not r.unpack, that is a completely different workflow that has nothing to do with reprojection.

Markus M

Unpack the DEM in the current location/mapset;
r.unpack n45w121.pack
WARNING: Difference between PROJ_INFO file of packed map and of current
location:
name: NAD83(HARN) / Oregon North
datum: nad83harn
ellps: grs80
proj: lcc
lat_1: 46
lat_2: 44.33333333333334
lat_0: 43.66666666666666
lon_0: -120.5
x_0: 2500000
y_0: 0
no_defs: defined

  • towgs84: 0.000,0.000,0.000
    ERROR: Projection of dataset does not appear to match current location. In
    case of no significant differences in the projection definitions,
    use the -o flag to ignore them and use current location definition.

If I correctly understand gdalwarp it functions like r.proj in reprojecting
a raster map from one location to another location. Here, both the source
and target locations have the same projection. Will gdalwarp work in this
case? Isn’t r.pack/r.unpack the correct module pair?

Please advise me of what I should do now or what other information I can
provide.

Regards,

Rich

grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

11

On Sun, 15 Sep 2019, Markus Metz wrote:

You should now use
r.proj loc=Oregon map=topography in=n45w121
and not r.unpack, that is a completely different workflow that has nothing
to do with reprojection.

Markus M,

I've tried two approaches to getting the 10m DEM (upper, left corner at 46N,
123W) into the project location: directly reprojecting from the import
location (in lon-lat decimal degrees), and unpacking from the r.packed map
in the statewide location (which successfully was projected from the import
location).

The locations are apparently not the problem; the statewide and project
locations are the same. I changed the project region to the state boundary
(the same region as the statewide location) and this makes no difference.

Both approaches fail; the resulting map has only 500 columns while the
source map has 11k+ columns.

Stepping back from this problem I recognized that the furthest western
longitude for the DEM immediately south of this problem one is too far to
the east at 121W. So, I can't use the 10m DEMs for lack of coverage.

I'm going back to the LiDAR data (which covers the entire north coast) and
can always resample to 5m or 10m if 1m is too fine for analyzing sediment
transport dynamics and fish populations in the entire river network.

I apppreciate your comments as well as those by Helmut and Markus N.

Best regards,

Rich