180 meridian and r.proj

[awickert]

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...?

[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.

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

[hamish]

Replying to 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

[awickert]

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.

[awickert]

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

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

[awickert]

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

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)