[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.
> >>
> >> 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.
> >>
> >> 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
> >>>>
> >>>>
> >>>>
> >>>>
> >>
> >>
> >
> >
> >
>
>
--
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
More information about the grass-user
mailing list