[GRASS-user] r.proj: Regions Do Not Match

[Hamish]

Rich wrote:
>   I am trying to reproject all
> data to a single projection: Albers Equal
> Area (NAD83) and g.region is:
> 
> projection: 99 (Albers Equal Area)
> zone:       0
> datum:      nad83
> ellipsoid:  grs80
> north:      43
> south:      34
> west:       -121
> east:       -112
> nsres:      1
> ewres:      1
> rows:       9
> cols:       9
> cells:      81
> 
>   The DEMs were originally in lat/lon (NAD83) and
> imported to that location.

if the originals are lat/lon then they need to be imported to a lat/lon
location, not a projected Albers location. pull them in from the lat/lon
location to the AEA location with r.proj if needed. I think you've already
done this, but the region bounds above don't look very Albers-like.

Albers will be in meters (or feet) -- Cartesian units.
Lat/lon will be in degrees -- angular units.
those don't mix without a projection.


> GRASS 6.5.svn (Nevada-aea):~/grassdata > r.proj in=dem_area15
> location=Nevada-ll mapset=PERMANENT method=cubic_f
...
> ERROR: Input raster map is outside current region
 
>   Running r.proj -p produces:
> 
> GRASS 6.5.svn (Nevada-aea):~/grassdata > r.proj -p in=dem_area15
> location=Nevada-ll mapset=PERMANENT method=cubic_f
> Input Projection Parameters:  +proj=longlat +no_defs
> +a=6378137
> +rf=298.257222101 +towgs84=0,0,0,0,0,0,0
> Input Unit Factor: 1
> Output Projection Parameters:  +proj=aea +lat_1=29.5
> +lat_2=45.5 +lat_0=23
> +lon_0=-96 +x_0=0 +y_0=0 +no_defs +a=6378137
> +rf=298.257222101
> +towgs84=0,0,0,0,0,0,0
> Output Unit Factor: 1
> Input map <dem_area15 at PERMANENT> in location
> <Nevada-ll>:
> Source cols: 3600
> Source rows: 3600
> Local north: 2064411.8954923
> Local south: 1973182.43184852
> Local west: -1787211.15607592
> Local east: -1679987.035036
> 
>   Is this a result of user error? If so, do I
> re-import these data? If not,
> what do I need to do to get all these data into a single
> location?

ok, that looks good, although I wonder why the east and west have become
negative.

r.proj pulls into the current region bounds, so try this step first, taking
bounds from "r.proj -p" output above: ('r.proj -g' is handy for that)

g.region n=2064411.8954923 s=1973182.43184852 e=-1679987.035036 \
   -p w=-1787211.15607592 rows=3600 cols=3600

which gives you:
nsres:      25.34151768
ewres:      29.78447807

maybe then try to round the resolution and bounds:
  g.region -a -p nsres=25 ewres=30

or both to 30m (or ft) if you prefer nice squares.


usually I try to have no less than the same number of rows and columns as
the original map, and if the convergence angle between the two projections
is large (g.region -n) add more rows+cols to compensate for the rotation.


once you are happy with the target resolution/bounds, only then run
r.proj to pull in the other map -into the current region box-.


you're very close,
Hamish