[GRASSLIST:2689] stitching geotiffs

GRASS GRASS user
,
1

Hi,

I am a newbie to GRASS (XT, Cygwin/XFree86; GRASS5.0) and GIS trying to stitch together a set of NLCD Geo-Tiff (USGS Landsat-5, etc.) data files for a set of contiguous states in the US. I finally figured out how to import them using r.in.gdal -o. (If there is some inherent problem with overriding the projection checking, please let me know.)

For starters, how do I go about importing the individual files with their corresponding UTM zone designations into the same project?

Assuming I can figure that out, do I use r.mask to remove the background from each of these maps?

Thanks for your help.

Martin
mdusaire@umn.edu

2

Martin du Saire wrote:

Hi,

I am a newbie to GRASS (XT, Cygwin/XFree86; GRASS5.0) and GIS trying to stitch together a set of NLCD Geo-Tiff (USGS Landsat-5, etc.) data files for a set of contiguous states in the US. I finally figured out how to import them using r.in.gdal -o. (If there is some inherent problem with overriding the projection checking, please let me know.)

If you are importing an image in UTM Zone X into a location defined as UTM Zone X, then you are doing the correct thing by using -o. (You should not use -o to import Zone X into a location of Zone Y). You may also want to use -e, which extends the location if needed.

For starters, how do I go about importing the individual files with their corresponding UTM zone designations into the same project?

Do I understand correctly that you have images in Zone X, Y Z ... and want to have them all in a single projection? If I understand correctly, then you need to use r.proj to reproject from the individual zones into the single target projection. But maybe I do not understand you goal.

Assuming I can figure that out, do I use r.mask to remove the background from each of these maps?

You probably do not need to mask. I am guessing that the background is NULLs?

--
Richard Greenwood
www.greenwoodmap.com

3

At 04:01 PM 2/19/2004 -0700, you wrote:

Martin du Saire wrote:

Hi,
I am a newbie to GRASS (XT, Cygwin/XFree86; GRASS5.0) and GIS trying to stitch together a set of NLCD Geo-Tiff (USGS Landsat-5, etc.) data files for a set of contiguous states in the US. I finally figured out how to import them using r.in.gdal -o. (If there is some inherent problem with overriding the projection checking, please let me know.)

If you are importing an image in UTM Zone X into a location defined as UTM Zone X, then you are doing the correct thing by using -o. (You should not use -o to import Zone X into a location of Zone Y). You may also want to use -e, which extends the location if needed.

So I can import a Zone X map into a Zone Y location? Is this any better/worse than using r.proj? The bulk of my area of interest is in a single UTM zone (Zone Y), so maybe the distortions generated by importing Zone X and Zone Z maps into a location of Zone Y will be small??

For starters, how do I go about importing the individual files with their corresponding UTM zone designations into the same project?

Do I understand correctly that you have images in Zone X, Y Z ... and want to have them all in a single projection? If I understand correctly, then you need to use r.proj to reproject from the individual zones into the single target projection. But maybe I do not understand you goal.

Yes, that is what I am attempting to do. I am stitching together maps for MN, IA (Zone 15), WI (Zone 16), and ND, SD (Zone 14). About 80% of my area of interest is in Zone 15.
I also came across a thread from Feb97 that suggested making symbolic links of mapsets and then running g.mapsets. What you suggest sound more straightforward. Would these procedures be equivalent?

Assuming I can figure that out, do I use r.mask to remove the background from each of these maps?

You probably do not need to mask. I am guessing that the background is NULLs?

The background is set to "0". Would a NULL be transparent? Maybe use: r.null map="name" setnull=0

--
Richard Greenwood
www.greenwoodmap.com

Thanks for your help.

Martin

4

Martin du Saire wrote:

> > I am a newbie to GRASS (XT, Cygwin/XFree86; GRASS5.0) and GIS trying to
> > stitch together a set of NLCD Geo-Tiff (USGS Landsat-5, etc.) data files
> > for a set of contiguous states in the US. I finally figured out how to
> > import them using r.in.gdal -o. (If there is some inherent problem with
> > overriding the projection checking, please let me know.)
>
> If you are importing an image in UTM Zone X into a location defined as UTM
> Zone X, then you are doing the correct thing by using -o. (You should not
> use -o to import Zone X into a location of Zone Y). You may also want to
> use -e, which extends the location if needed.

So I can import a Zone X map into a Zone Y location?

You can (i.e. nothing will actually stop you), but you shouldn't. You
may have misread or misinterpreted what Richard wrote.

Is this any better/worse than using r.proj?

Worse.

The bulk of my area of interest is in a
single UTM zone (Zone Y), so maybe the distortions generated by importing
Zone X and Zone Z maps into a location of Zone Y will be small??

That depends upon whether you consider an error of a few tens of
kilometres "small".

>>For starters, how do I go about importing the individual files with their
>>corresponding UTM zone designations into the same project?
>
>Do I understand correctly that you have images in Zone X, Y Z ... and want
>to have them all in a single projection? If I understand correctly, then
>you need to use r.proj to reproject from the individual zones into the
>single target projection. But maybe I do not understand you goal.

Yes, that is what I am attempting to do. I am stitching together maps for
MN, IA (Zone 15), WI (Zone 16), and ND, SD (Zone 14). About 80% of my area
of interest is in Zone 15.

In which case, you need three locations, one for each zone. Import
each map into the correct location, then start GRASS in the zone 15
location and use r.proj to re-project the zone 14/16 maps.

I also came across a thread from Feb97 that suggested making symbolic links
of mapsets and then running g.mapsets. What you suggest sound more
straightforward. Would these procedures be equivalent?

No. If the maps use different projections, you have to use r.proj at
some point.

>>Assuming I can figure that out, do I use r.mask to remove the background
>>from each of these maps?
>
>You probably do not need to mask. I am guessing that the background is NULLs?

The background is set to "0". Would a NULL be transparent?

It depends upon context; for r.patch, NULL is effectively
"transparent", i.e. non-NULL cells overwrite whatever is "beneath"
them, NULL cells don't.

Maybe use: r.null map="name" setnull=0

Yes.

--
Glynn Clements <glynn.clements@virgin.net>

5

Martin du Saire wrote:

At 04:01 PM 2/19/2004 -0700, you wrote:

Martin du Saire wrote:

Hi,
I am a newbie to GRASS (XT, Cygwin/XFree86; GRASS5.0) and GIS trying to stitch together a set of NLCD Geo-Tiff (USGS Landsat-5, etc.) data files for a set of contiguous states in the US. I finally figured out how to import them using r.in.gdal -o. (If there is some inherent problem with overriding the projection checking, please let me know.)

If you are importing an image in UTM Zone X into a location defined as UTM Zone X, then you are doing the correct thing by using -o. (You should not use -o to import Zone X into a location of Zone Y). You may also want to use -e, which extends the location if needed.

So I can import a Zone X map into a Zone Y location? Is this any better/worse than using r.proj? The bulk of my area of interest is in a single UTM zone (Zone Y), so maybe the distortions generated by importing Zone X and Zone Z maps into a location of Zone Y will be small??

No. Sorry I was not clear. You must use r.proj. Import Zone 14 maps into a Zone 14 location. Zone 15 into a Zone 15 location, etc. And then use r.proj to bring your Zone 14, and Zone 16 maps into your Zone 15 location.

For starters, how do I go about importing the individual files with their corresponding UTM zone designations into the same project?

Do I understand correctly that you have images in Zone X, Y Z ... and want to have them all in a single projection? If I understand correctly, then you need to use r.proj to reproject from the individual zones into the single target projection. But maybe I do not understand you goal.

Yes, that is what I am attempting to do. I am stitching together maps for MN, IA (Zone 15), WI (Zone 16), and ND, SD (Zone 14). About 80% of my area of interest is in Zone 15.
I also came across a thread from Feb97 that suggested making symbolic links of mapsets and then running g.mapsets. What you suggest sound more straightforward. Would these procedures be equivalent?

I am not familiar with g.mapsets so I can not comment.

Assuming I can figure that out, do I use r.mask to remove the background from each of these maps?

You probably do not need to mask. I am guessing that the background is NULLs?

The background is set to "0". Would a NULL be transparent? Maybe use: r.null map="name" setnull=0

Yes, that sound correct.

What do you want to get as a final result? A single map (image) or a whole bunch of maps (image tiles)? If you want a single map, then you would probably want to use r.patch to assemble all the individual maps in the Zone 15 location.

Rich

--
Richard Greenwood
www.greenwoodmap.com

6

Outstanding! Thank you Richard and Glynn. I think I may eventually have to r.patch the individual maps. I was playing around with two maps within a single zone, and the r.stats -a -n map1,map2 output did not reflect what I expected to be the composite of the two (Actually, the numbers made no sense at all). One last question:

Without having specific datum in my adjacent zones (14, 16), is there some way to make upper-bound estimates on the amount of distortion resulting from the r.proj method=nearest (the maps are land-use/cover)?

Thanks again.

Martin

At 10:06 AM 2/20/2004 -0700, you wrote:

Martin du Saire wrote:

At 04:01 PM 2/19/2004 -0700, you wrote:

Martin du Saire wrote:

Hi,
I am a newbie to GRASS (XT, Cygwin/XFree86; GRASS5.0) and GIS trying to stitch together a set of NLCD Geo-Tiff (USGS Landsat-5, etc.) data files for a set of contiguous states in the US. I finally figured out how to import them using r.in.gdal -o. (If there is some inherent problem with overriding the projection checking, please let me know.)

If you are importing an image in UTM Zone X into a location defined as UTM Zone X, then you are doing the correct thing by using -o. (You should not use -o to import Zone X into a location of Zone Y). You may also want to use -e, which extends the location if needed.

So I can import a Zone X map into a Zone Y location? Is this any better/worse than using r.proj? The bulk of my area of interest is in a single UTM zone (Zone Y), so maybe the distortions generated by importing Zone X and Zone Z maps into a location of Zone Y will be small??

No. Sorry I was not clear. You must use r.proj. Import Zone 14 maps into a Zone 14 location. Zone 15 into a Zone 15 location, etc. And then use r.proj to bring your Zone 14, and Zone 16 maps into your Zone 15 location.

<snip>

7

Martin du Saire wrote:

Outstanding! Thank you Richard and Glynn. I think I may eventually have to r.patch the individual maps. I was playing around with two maps within a single zone, and the r.stats -a -n map1,map2 output did not reflect what I expected to be the composite of the two (Actually, the numbers made no sense at all). One last question:

Without having specific datum in my adjacent zones (14, 16), is there some way to make upper-bound estimates on the amount of distortion resulting from the r.proj method=nearest (the maps are land-use/cover)?

Thanks again.

Martin

That's a hard question that I will not try to answer directly. Map distortion can involve linear, directional, and/or area distortion. So first you have to define which of those are relevant to your application. I think in your case the most significant distortion will be angular. And I would guess that areas would be most important for land use applications. Within a UTM zone I am thinking that area distortion is limited to about 0.04% because the central meridian has a scale factor of 0.9996 (anyone is welcome to correct me on that). But as you move outside the zone, the increase in distortion would not be linear.

I am kind of a trial-and-error guy, so if it was me, I would simply reproject the extreme cases into the target location and see if they meet your requirements. Extending UTM zones is not generally considered a good practice (although I am guilty of having done it a couple times). A better solution might be to create a projection for your project.

--
Richard Greenwood
www.greenwoodmap.com

8

On Fri, 20 Feb 2004, Richard Greenwood wrote:

Martin du Saire wrote:

> Outstanding! Thank you Richard and Glynn. I think I may eventually
> have to r.patch the individual maps. I was playing around with two maps
> within a single zone, and the r.stats -a -n map1,map2 output did not
> reflect what I expected to be the composite of the two (Actually, the
> numbers made no sense at all). One last question:
>
> Without having specific datum in my adjacent zones (14, 16), is there
> some way to make upper-bound estimates on the amount of distortion
> resulting from the r.proj method=nearest (the maps are land-use/cover)?
>
> Thanks again.
>
> Martin

That's a hard question that I will not try to answer directly. Map
distortion can involve linear, directional, and/or area distortion. So
first you have to define which of those are relevant to your
application. I think in your case the most significant distortion will
be angular. And I would guess that areas would be most important for
land use applications. Within a UTM zone I am thinking that area
distortion is limited to about 0.04% because the central meridian has a
scale factor of 0.9996 (anyone is welcome to correct me on that). But as
you move outside the zone, the increase in distortion would not be linear.

I am kind of a trial-and-error guy, so if it was me, I would simply
reproject the extreme cases into the target location and see if they
meet your requirements. Extending UTM zones is not generally considered
a good practice (although I am guilty of having done it a couple times).
A better solution might be to create a projection for your project.

For visual checking you can use r.mapcalc to create a raster grid
in one location, then project it into an other location and check
the distorion. e.g.

  r.mapcalc testgrid='if((row()%5) && (col()%5))'

will create a grid where every 5th row and column have a different color.
I bet someone has a more sofisticated recipe.

9

At 03:38 PM 2/20/2004 -0700, you wrote:

Martin du Saire wrote:

Outstanding! Thank you Richard and Glynn. I think I may eventually have to r.patch the individual maps. I was playing around with two maps within a single zone, and the r.stats -a -n map1,map2 output did not reflect what I expected to be the composite of the two (Actually, the numbers made no sense at all). One last question:
Without having specific datum in my adjacent zones (14, 16), is there some way to make upper-bound estimates on the amount of distortion resulting from the r.proj method=nearest (the maps are land-use/cover)?
Thanks again.
Martin

That’s a hard question that I will not try to answer directly. Map distortion can involve linear, directional, and/or area distortion. So first you

Thanks again.

Having a little trouble with r.proj though:

ERROR> Input map is outside current region

I read the thread that involved some problem with southern hemisphere projection data from U of Maryland…I am working in the NE USA…geotiffs from the USGS site (for what that is worth).

Region/projection info for the location I am attempting to project the map into are:

g.region -p

projection: 1 (UTM)
zone: 19
datum: nad83
ellipsoid: grs80
north: 2734830
south: 1944930
west: 1721790
east: 2041380
nsres: 30
ewres: 30
rows: 26330
cols: 10653

g.projinfo

PROJ_INFO file:
name: UTM
datum: nad83
dx: 0.000000
dy: 0.000000
dz: 0.000000
proj: utm
ellps: grs80
a: 6378137.0000000000
es: 0.0066943800
f: 298.2572221010
zone: 19

PROJ_UNITS file:
unit: meter
units: meters
meters: 1.0

Region/projection for the map are:

projection: 1 (UTM)
zone: 18
datum: nad83
ellipsoid: grs80
north: 2368230
south: 1944930
west: 1721790
east: 1988670
nsres: 30
ewres: 30
rows: 14110
cols: 8896

PROJ_INFO file:
name: UTM
datum: nad83
dx: 0.000000
dy: 0.000000
dz: 0.000000
proj: utm
ellps: grs80
a: 6378137.0000000000
es: 0.0066943800
f: 298.2572221010
zone: 18

PROJ_UNITS file:
unit: meter
units: meters
meters: 1.0

I tried to the expand the current region a little (North:280000, South:1900000, West:1700000, East:2100000), but that didn’t help.

Help. Thanks again, this GIS stuff is quite fun.

Martin

10

Martin du Saire wrote:

At 03:38 PM 2/20/2004 -0700, you wrote:

Martin du Saire wrote:

Outstanding! Thank you Richard and Glynn. I think I may eventually have to r.patch the individual maps. I was playing around with two maps within a single zone, and the r.stats -a -n map1,map2 output did not reflect what I expected to be the composite of the two (Actually, the numbers made no sense at all). One last question:
Without having specific datum in my adjacent zones (14, 16), is there some way to make upper-bound estimates on the amount of distortion resulting from the r.proj method=nearest (the maps are land-use/cover)?
Thanks again.
Martin

That's a hard question that I will not try to answer directly. Map distortion can involve linear, directional, and/or area distortion. So first you

Thanks again.

Having a little trouble with r.proj though:

ERROR> Input map is outside current region

You are on the right track - enlarge your region to encompass the area of the maps that you are bringing in. And then set your current region to your default region (g.region 1). Remember that GRASS only works within the bounds of your current region, is if you had zoomed into an area, that's the only area in the target location, that would be the only area that would be processed with most GRASS commands.

Rich

--
Richard Greenwood
www.greenwoodmap.com

11

Martin du Saire wrote:

>>Outstanding! Thank you Richard and Glynn. I think I may eventually
>>have to r.patch the individual maps. I was playing around with two maps
>>within a single zone, and the r.stats -a -n map1,map2 output did not
>>reflect what I expected to be the composite of the two (Actually, the
>>numbers made no sense at all). One last question:
>>Without having specific datum in my adjacent zones (14, 16), is there
>>some way to make upper-bound estimates on the amount of distortion
>>resulting from the r.proj method=nearest (the maps are land-use/cover)?
>>Thanks again.
>>Martin
>
>That's a hard question that I will not try to answer directly. Map
>distortion can involve linear, directional, and/or area distortion. So
>first you

Thanks again.

Having a little trouble with r.proj though:

ERROR> Input map is outside current region

I read the thread that involved some problem with southern hemisphere
projection data from U of Maryland...I am working in the NE USA...geotiffs
from the USGS site (for what that is worth).

Region/projection info for the location I am attempting to project the map
into are:

g.region -p

projection: 1 (UTM)
zone: 19
datum: nad83
ellipsoid: grs80
north: 2734830
south: 1944930
west: 1721790
east: 2041380

That's way off. Valid eastings for UTM are within 350km of the false
500km false easting, i.e. ~166000-834000. That's at the equator; the
valid range decreases with increasing latitude. So, your region is off
by a couple of zones.

Zone 18 is 72-78 degrees West, while zone 19 is 66-72 degrees West.

Also, that particular region appears way South (somewhere to the East
of Guadeloupe); The North value of 2734830 is about right for the
southern tip of Florida (but that would be zone 17).

--
Glynn Clements <glynn.clements@virgin.net>