[GRASS-dev] [GRASS GIS] #1594: 180 meridian and r.proj

GRASS GIS GRASS GIS developer
1

#1594: 180 meridian and r.proj
------------------------------+---------------------------------------------
Reporter: awickert | Owner: grass-dev@…
     Type: defect | Status: new
Priority: normal | Milestone:
Component: Default | Version: 6.4.1
Keywords: r.proj, meridian | Platform: Linux
      Cpu: x86-64 |
------------------------------+---------------------------------------------
On versions 6.4.1. and 7.0, when I project a map that is in lat/lon
coordinates, only the region to the west of the 180 meridian (i.e. that up
to 180 degrees E) is projected and brought into the new coordinate system.
There must be some way of making r.proj realize that the rest of the map
exists...?

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/1594&gt;
GRASS GIS <http://grass.osgeo.org>

2

#1594: 180 meridian and r.proj
------------------------------+---------------------------------------------
Reporter: awickert | Owner: grass-dev@…
     Type: defect | Status: new
Priority: normal | Milestone:
Component: Raster | Version: 6.4.1
Keywords: r.proj, meridian | Platform: Linux
      Cpu: x86-64 |
------------------------------+---------------------------------------------
Changes (by hamish):

  * component: Default => Raster

Comment:

hmmm, I've never seen this problem with raster maps and I'm at 170E
longitude. If a reprojection hits a singularity typically you know about
it immediately with a flood of warning messages.

can you provide more details on how to reproduce the problem and possibly
some test data?

what projection system are you reprojecting into?

do you know if the original dataset was 0-360 or -180 to +180?

thanks,
Hamish

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1594#comment:1&gt;
GRASS GIS <http://grass.osgeo.org>

3

#1594: 180 meridian and r.proj
------------------------------+---------------------------------------------
Reporter: awickert | Owner: grass-dev@…
     Type: defect | Status: new
Priority: normal | Milestone:
Component: Raster | Version: 6.4.1
Keywords: r.proj, meridian | Platform: Linux
      Cpu: x86-64 |
------------------------------+---------------------------------------------

Comment(by hamish):

Replying to [comment:1 hamish]:
> hmmm, I've never seen this problem with raster maps and I'm at
> 170E longitude. If a reprojection hits a singularity typically
> you know about it immediately with a flood of warning messages.

... which is typically only a problem with a one-way projection like
Winkel Tripel. Unlike r.proj, v.proj uses a forward projection and in
those cases the work-around is to do `r.to.vect feature=point`, v.proj,
then `v.out.ascii | r.in.xyz`. Can't recall if I ever wrote an addon
script to automate that.

Hamish

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1594#comment:2&gt;
GRASS GIS <http://grass.osgeo.org>

4

#1594: 180 meridian and r.proj
------------------------------+---------------------------------------------
Reporter: awickert | Owner: grass-dev@…
     Type: defect | Status: new
Priority: normal | Milestone:
Component: Raster | Version: svn-trunk
Keywords: r.proj, meridian | Platform: Linux
      Cpu: x86-64 |
------------------------------+---------------------------------------------
Changes (by awickert):

  * version: 6.4.1 => svn-trunk

Comment:

Hi Hamish,

Sorry that I completely overlooked this bug report and let it fall by the
wayside!

Here is exactly what I am doing:

1. Importing a global gridded data set into a GRASS GIS location (we'll
call it "L1"; this is unprojected: WGS84 = EPSG 4326)

2. Creating a new GRASS GIS location (L2), in the same unprojected
coordinate system (WGS84)

3. Inside L2, setting the region to extend across the 180 meridian. In my
case, here is the computational region, which is at the same resolution as
the global gridded data:

> g.region -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 80N
south: 50N
west: 155E
east: 125W
nsres: 0:00:30
ewres: 0:00:30
rows: 3600
cols: 9600
cells: 34560000

4. I run the following command to import the data set, in this case, the
GEBCO_08 30-arcsecond global digital elevation model:

> r.proj loc=GlobalDatabase30as in=toponow

5. I look at the output data extent - only the part of the region to the
west of the 180 meridian is imported.

I could upload the raster, but it is huge... and this is reproducible by
just creating a map of zeros with r.mapcalc in the global location ("L1")
and trying to bring it into L2 using r.proj.

Hope this helps - sorry about the *long* delay - and happy to help with
some coding if that is required to take care of this!

Andy

PS - I just uploaded to today's svn-trunk, so the problem is current.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/1594#comment:3&gt;
GRASS GIS <http://grass.osgeo.org>

5

#1594: 180 meridian and r.proj
------------------------------+---------------------------------------------
Reporter: awickert | Owner: grass-dev@…
     Type: defect | Status: new
Priority: normal | Milestone:
Component: Raster | Version: svn-trunk
Keywords: r.proj, meridian | Platform: Linux
      Cpu: x86-64 |
------------------------------+---------------------------------------------

Comment(by awickert):

Here is my current workaround, using r.proj and r.patch:

{{{
#!python

from grass import script as grass
import re

inloc = 'GlobalDatabase30as'
non_decimal = re.compile(r'[^\d.]+') # for my hack-ish solution below
# "toponow" NEEDS to be added (some other good too):
basicfiles = ['precip_2000_2004', 'et_2000_2004', 'cellsize_km2',
'cellsize_meters2', 'zeros', 'toponow']
for infile in basicfiles:
   print infile
   # Here we make sure we are using the same resolution as the input map.
   # Getting the resolution - there has got to be a much better way to do
this,
   # but this is the first thing that came to mind
   shreg = grass.parse_command('r.proj', location=inloc,
input='precip_2000_2004', output='tmp0', overwrite=True,
flags='g')['n'].split(' ')[-2:]
   for i in range(len(shreg)):
     shreg[i] = int(non_decimal.sub('', (shreg[i])))
   grass.run_command('g.region', nsres=shreg[0]/180., ewres=shreg[1]/360.,
e='180E')
   grass.run_command('r.proj', location=inloc, input=infile, output='tmp0',
overwrite=True)
   grass.run_command('g.region', region='default')
   grass.run_command('g.region', sres=shreg[0]/180., ewres=shreg[1]/360.,
w='180W')
   grass.run_command('r.proj', location=inloc, input=infile, output='tmp1',
overwrite=True)
   grass.run_command('g.region', region='default')
   grass.run_command('r.patch', input='tmp0,tmp1', output=infile,
overwrite=True)
}}}

Cheers,

Andy Wickert

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/1594#comment:4&gt;
GRASS GIS <http://grass.osgeo.org>

6

#1594: 180 meridian and r.proj
------------------------------+---------------------------------------------
Reporter: awickert | Owner: grass-dev@…
     Type: defect | Status: new
Priority: normal | Milestone:
Component: Raster | Version: svn-trunk
Keywords: r.proj, meridian | Platform: Linux
      Cpu: x86-64 |
------------------------------+---------------------------------------------

Comment(by awickert):

A few mistakes in my previously-posted code; use this instead:
Cheers,
Andy

{{{
#!python

from grass import script as grass
import re

grass.run_command('g.region', save='default')

inloc = 'GlobalDatabase30as'
non_decimal = re.compile(r'[^\d.]+') # for my hack-ish solution below
# "toponow" NEEDS to be added (some other good too):
basicfiles = ['precip_2000_2004', 'et_2000_2004', 'cellsize_km2',
'cellsize_meters2', 'zeros', 'toponow']
for infile in basicfiles:
   print infile
   # Here we make sure we are using the same resolution as the input map.
   # Getting the resolution - there has got to be a much better way to do
this,
   # but this is the first thing that came to mind
   shreg = grass.parse_command('r.proj', location=inloc, input=infile,
output='tmp0', overwrite=True, flags='g')['n'].split(' ')[-2:]
   for i in range(len(shreg)):
     shreg[i] = int(non_decimal.sub('', (shreg[i])))
   grass.run_command('g.region', nsres=180./shreg[0], ewres=360./shreg[1],
e='180E')
   grass.run_command('r.proj', location=inloc, input=infile, output='tmp0',
overwrite=True)
   grass.run_command('g.region', region='default')
   grass.run_command('g.region', nsres=180./shreg[0], ewres=360./shreg[1],
w='180W')
   grass.run_command('r.proj', location=inloc, input=infile, output='tmp1',
overwrite=True)
   grass.run_command('g.region', region='default')
   grass.run_command('r.patch', input='tmp0,tmp1', output=infile,
overwrite=True)

}}}

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/1594#comment:5&gt;
GRASS GIS <http://grass.osgeo.org>