[GRASSLIST:7996] Re: r.proj takes a very long time / does not complete
David Finlayson
david.p.finlayson at gmail.com
Fri Aug 19 19:03:01 EDT 2005
The only reason to project the corners is to allow you to make a guess
at an appropriate resolution for your data before you run r.proj.
Another way would be to run r.proj on the data and deliberately use a
low-resolution setting (just to get the raster into projected space
quickly). You can then use d.measure to measure the length and width
of the data area (in meters or feet now instead of fractions of a
degree). Then run r.proj again at a better resolution for your data.
I'll get on my soap box...
Take a few examples, the National DEM of the United States (USGS) is
distributed in "geographic" projection, which is ESRI-speak for it has
no projection at all.
At 30N the length of 1 arc second is about 27 m, at 50N it is 20 m. So
the longitudinal "resolution" of the US Map narrows by about 25% from
the bottom of the US to the top. Why they chose to do this is beyond
me. But at least in the figures of the USA they show a beautiful
Lambert projection of the US. Clearly, THEY know you need to project
the data before you can use it. (http://ned.usgs.gov/)
The national coastal relief model does the same thing (storing data in
lat-long grids) with their bathymetric data. Only here, the maps they
produced used a Mercator projection (the default projection of all
oceanographers) which grossly distorts the shape of North America:
http://www.ngdc.noaa.gov/mgg/coastal/coastal.html
My favorite map is this one which manages to distort the longitudinal
distance of the Pacific Northwest by 156%:
http://www.ngdc.noaa.gov/mgg/coastal/grddas08/grddas08.htm
With cartography like this, what else lurks in these data?
OK, I'm standing down.
David
On 8/19/05, Jason Horn <jhorn at bu.edu> wrote:
> David,
>
> Thanks for your thoughtful comments. I think I understand what you
> are saying about resolution not really existing in a lat-lon
> coordinate system. I'll try your idea of projecting the 4 corners,
> and then computing the resolution. One question, isn't this what
> r.proj would do anyway?
>
> I would like to see your spreadsheet. If you can, please email it to
> jhorn at bu.edu.
>
> Thanks
>
> - Jason
>
>
> On Aug 19, 2005, at 4:44 AM, David Finlayson wrote:
>
> > 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
> >
> >
> >
>
>
--
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