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