[GRASSLIST:7970] Re: r.proj takes a very long time / does not complete

Morten Hulden morten at untamo.net
Thu Aug 18 13:14:54 EDT 2005

Jason Horn wrote:

> Thanks so much for pointing out the problem.  I have seceded in  
> projecting the raster layer, with the help of the v.in.region box  
> trick.  I have set the resolution to an approximation of what I think  
> is right.  However, how can I really match the resolution in the new  
> state plane location to the resolution from the old lat-lon location?

Since latlon resolution is measured in degrees while most projections 
use meters unit you have to approximate. At the equator one degree = 
40,000,000/360 meters for both longitude and latitude but further north 
or south the length of a one degree longitude arc decreases since 
meridians converge toward poles.

But assuming a fixed value 40,000,000/360 for a one degree arc then:
the resolution of your latlon map 0:00:38.190402 becomes:

40,000,000 / 360 / 60 / 60 * 38.190402 = 1178.7161111111 meters

You may of course use a more exact value for the earth's perimeter 
instad of 40,000,000 meters but since the projection involves resampling 
  errors anyway I don't think it will make a big difference. 1000 meters 
resolution would probably be close enough.

> On Aug 15, 2005, at 5:25 PM, Ian MacMillan wrote:
>> A simple way to do this is to make a vector box of your region in  the 
>> first location using v.in.region.  Go to your new location,  v.proj 
>> that vector in, and set your region to the new vector.  Then  use 
>> r.proj on the raster you want to project.

>> On Aug 15, 2005, at 2:14 PM, Morten Hulden wrote:
>>> You should adjust the boundaries or resolution of your NAD83  
>>> location to something more reasonable. Before running r.proj try  to 
>>> set a resolution value that roughly matches the resolution of  your 
>>> input map (g.region res= ), e.g. something like 1000 m.
>>> Don't know about boundaries for your location (TX unknown area to  
>>> me) but you may want to check those as well. No need to cover more  
>>> than what is necessary.

I don't think this is necessary to use the v.in.region trick here 
because r.proj will trim down the map anyway to the area that is 
overlapping the current region. What I meant was that if you have a very 
large region in your location (847325 x 622527 cells) you probably do 
not want to have a small resolution of 1 meter. Either reduce the region 
to an area of interest or increase the resolution, or both.

My own preferrence is to use a rather coarse resolution for the default 
region within a location, say a utm zone within my country - 456 meters 
matches both landsat mosaics and 30-arcsec DEMs close enough. For 
selected areas of interest within the location I create regions with 
decreased coverage but increased resolution, 57, 28.5 (landsat5), and 
14.25 meters (landsat7 visual band) and even less than that for my own 
house and back garden.

For maps with greater resolution than my default region's I create a 
region with matching resolution before projecting. I gain nothing by 
trying to project a coarse map into a finer grained region, I only lose 
disk space and time.

This way I can keep all maps I want to analyze together within a single 
location in Grass without wasting space and computation time on areas 
not relevant for the tasks I want to do.


>>> Jason Horn wrote:
>>>> Morten,
>>>> Thanks for looking at this for me.  Below are the files you asked  for.
>>>> PROJ_INFO for the source lat-lon location:
>>>> name: Latitude-Longitude
>>>> datum: wgs84
>>>> towgs84: 0.000,0.000,0.000
>>>> proj: ll
>>>> ellps: wgs84
>>>> cellhd file for  KEWX20040710_025420:
>>>> proj:       3
>>>> zone:       0
>>>> north:      34:00:17.832478N
>>>> south:      25:11:59.799175N
>>>> east:       92:43:25.599215W
>>>> west:       103:19:56.000785W
>>>> cols:       1000
>>>> rows:       830
>>>> e-w resol:  0:00:38.190402
>>>> n-s resol:  0:00:38.190402
>>>> format:     -1
>>>> compressed: 1
>>>> PROJ_INFO file for the destination NAD83 location:
>>>> name: Lambert Conformal Conic
>>>> datum: nad83
>>>> towgs84: 0.000,0.000,0.000
>>>> proj: lcc
>>>> ellps: grs80
>>>> a: 6378137.0000000000
>>>> es: 0.0066943800
>>>> f: 298.2572221010
>>>> lat_0: 31.1666670000
>>>> lat_1: 27.4166670000
>>>> lat_2: 34.9166670000
>>>> lon_0: -100.0000000000
>>>> x_0: 1000000.0000000000
>>>> y_0: 1000000.0000000000
>>>> WIND file for the destination NAD83 location:
>>>> proj:       99
>>>> zone:       0
>>>> north:      1196061
>>>> south:      573534
>>>> east:       1563147
>>>> west:       715822
>>>> cols:       847325
>>>> rows:       622527
>>>> e-w resol:  1
>>>> n-s resol:  1
>>>> top:        1
>>>> bottom:     0
>>>> cols3:      847325
>>>> rows3:      622527
>>>> depths:     1
>>>> e-w resol3: 1
>>>> n-s resol3: 1
>>>> t-b resol:  1
>>>> Thanks again!
>>>> - Jason
>>>> On Aug 13, 2005, at 12:14 PM, Morten Hulden wrote:
>>>>> Jason Horn wrote:
>>>>>> I have a long-standing problem with r.proj.  I have never been   
>>>>>> able to complete a successful re-projection of a raster file.   
>>>>>> For  example, I am trying to pull in a raster file from a  
>>>>>> location that  is lat-lon, WGS84 into a location that is state- 
>>>>>> plane meters,  NAD83.  The raster grid is relatively small,  1000 
>>>>>> cols by 830  rows.  I start r.proj like this:
>>>>>> r.proj input=KEWX20040710_025420 location=TX_ll_WGS84   
>>>>>> mapset=PERMANENT output=proj1 method=nearest Although my CPU  
>>>>>> stays  pegged at close to 100% usage, the progress bar never  
>>>>>> moves  forward and nothing happens, even if I leave it for many  
>>>>>> hours.   Can anyone tell me what's going on here?
>>>>> In order to help you I would like to see:
>>>>> 1) The PROJ_INFO file for your latlon location
>>>>> 2) The cellhd file of the KEWX20040710_025420 map in the latlon   
>>>>> location
>>>>> 3) The PROJ_INFO file for your NAD83 location
>>>>> 4) The WIND file for your NAD83 location

More information about the grass-user mailing list