Re: [GRASS-dev] [GRASS GIS] #1594: 180 meridian and r.proj
#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: 3456 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 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] [GRASS GIS] #1594: 180 meridian and r.proj
#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
GRASS GIS http://grass.osgeo.org
___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] [GRASS GIS] #1594: 180 meridian and r.proj
#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
GRASS GIS http://grass.osgeo.org
___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev
