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

David Finlayson david.p.finlayson at gmail.com
Fri Aug 19 04:44:48 EDT 2005

```I think that it is a conceptual mistake to consider a raster data set
in lat-long to have a "resolution" in the same sense as a projected
map. Degrees of latitude and longitude do not have constant spacing.
So a single resolution defined in meters (or whatnot) doesn't exist
(even on the globe) Now add to that, GIS can't work on ellipsoidal
surfaces and must project the globe onto a flat plane. So not only do
degrees have variable spacing to begin with, they must be warped to
get them into the flat-earth of the GIS.

I think it is best to consider lat-long grids a convenient storage
mechanism for variably spaced points because by the time you can
actually use them in a GIS (after projection and conversion to a
Cartesian coordinate system), that's what they will be. The point is
to remember that the conversion from geographic coordinates to
projected coordinates is a major transformation and not a simple
coordinate shift and resample problem.

As for how to estimate the resolution, I usually run the 4 corners of
my geographic data through proj to get the Cartesian coordinates
(State Plane, UTM, etc.). Now everything is on "flat-earth" and you
can calculate the appropriate resolution for your data with simple
arithmetic (after all, that's why we invented map projections!).

As a curiosity, if you want a spreadsheet that calculates the
ellipsoidal distance between any two points on the earth, send me an
email and I'll post it up on the web.

David

On 8/17/05, Jason Horn <jhorn at bu.edu> wrote:
> Ian, Morten,
>
> 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?
>
> - Jason
>
>
> 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.  You will still have
> > to set the resolution on your own though.  Hope this makes sense.
> >
> > -Ian
> >
> > On Aug 15, 2005, at 2:14 PM, Morten Hulden wrote:
> >
> >
> >>
> >>
> >> No wonder it takes forever. Your output location has 527482690275
> >> cells (847325 x 622527) to project.
> >>
> >> 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.
> >>
> >> 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
> >>> 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
> >>>>
> >>>>
> >>>>
> >>>>
> >>
> >>
> >
> >
> >
>
>

--
David Finlayson
Marine Geology & Geophysics
School of Oceanography
Box 357940
University of Washington
Seattle, WA  98195-7940
USA

Office: Marine Sciences Building, Room 112
Phone: (206) 616-9407
Web: http://students.washington.edu/dfinlays

```