Importing maps having different resolutions

rgp at abba rgp at abba
Thu Jul 14 12:44:30 EDT 1994


> From grass-lists-owner%max.cecer.army.mil at munnari.oz.au Thu Jul 14 08:02:48 1994
> Date: Wed, 13 Jul 1994 15:06:17 -0400 (EDT)
> From: "Robert Fraser (FO 1999)" <rfraser at minerva.cis.yale.edu>
> Subject: Importing maps having different resolutions
> Sender: grass-lists-owner at max.cecer.army.mil
> Reply-To: grassu-list at max.cecer.army.mil
> To: grassu-list at max.cecer.army.mil
> Mime-Version: 1.0
> Content-Type> : > TEXT/PLAIN> ; > charset=US-ASCII> 
> Content-Length: 751
> 
> I have successfully imported layers from IDRISI (soils, DEM, etc.) and 
> ERMAPPER (TM band) into GRASS using the r.in.ascii program.  All of these 
> layers have 30m resolution and the same number of rows and columns (i.e. 229 
> by 417).
> 
> I am now trying to import a Landsat MSS band having 56*79 m resolution 
> into the same LOCATION.  In ERMAPPER I windowed out an 87 row by 223 column 
> area corresponding to the previous layers and created an ascii file.  The 
> problem is, after I run r.in.ascii and try to display the band, I get a 
> series of discontinuous, coloured bands.
> 
> Does anyone know how I could resample my MSS band to the 30m resolution 
> and bring it into my 30m resolution LOCATION?
> 
> Thanks,
> 
> Robert Fraser
> Yale School of Forestry
> Hydro Lab
> 
It is unnecessary to convert to ascii to take a single band of a byte or integer
ERMapper image into Grass.You simply need to have each band as a separate image
and then run the shell script attached.
However because you are wanting to change the cell resolution also, the problem
becomes a little more complicated.
The following steps will achieve the required result.

Set up an algorithm to display the MSS band area of interest in ERMapper.

In ERMAPPER use the Algorithm Mode 'output to dataset' and select 'output to 
dataset options' in the main IMAGES menu.Click on 'Set Defaults'.

Examine the number of image lines and pixels which will be output by default.
If you require a square pixel recalculate the number of lines and pixels
for this area which would give the required cell size and enter them in
the 'Cells Across','Cells Down' fields.Select 'Output Data Type' to be
'Unsigned 8 bit Integer' or 'Signed 16 bit Integer' , Grass can handle either
format.The former is appropriate in this case.Select an 'Output Path Name'
for the new image and then hit 'GO' on the main menu.

The default resampling for the new cell size is nearest neighbour, OR if you
tick the 'Smooth Resampling' box in the algorithm window you will get bilinear
resampling.

Once you have produced your new image you are ready to take it into Grass.
To do this go into Grass select the appropriate MapSet and run the attached
shell script.
This script will build a Grass header file for your new image in cellhd 
by examining the ERMapper header and extracting the parts required.It then 
links the ERMapper image into your cell directory so that you don't have to 
have two copies on your disk.

#### Start Import Shell Script ####
-------------cut here
#!/bin/sh
#
# This shell script will create a GRASS cellhd file from an ER Mapper .ers
# file.  
# This shell script program shows
# how easy it is to copy from ER Mapper's ASCII format header files in to the
# GRASS header format.
#
if [ "$#" -ne "2" ]
then
echo "Use full path to ERMapper image from within a GRASS shell"
echo "Do not include the .ers extension"
echo "Usage: import_ermapper_to_grass <file> $LOCATION/cellhd/<file>"
exit 1
fi
ln -s $1 $LOCATION/cell/$2
ERS=$1.ers
GRASS=$LOCATION/cellhd/$2
echo > $GRASS "# GRASS $GRASS cellhd file created from ER Mapper $ERS file"
tr < $ERS '[a-z]' '[A-Z]' | awk -F= '
BEGIN {
        offsetx = 0.0;
        offsety = 0.0;
        xdim = 0.0;
        ydim = 0.0;
        regx = 0.0;
        regy = 0.0;
        nrows = 0;
        ncols = 0;
        bytorder = 0;
        form = 0;
}
/PROJECTION/               { proj = $2 }
/REGISTRATIONCELLX/        { offsetx = $2 }
/REGISTRATIONCELLY/        { offsety = $2 }
/METERSX/ || /EASTINGS/    { regx = $2 }
/METERSY/ || /NORTHINGS/   { regy = $2 }
/XDIMENSION/               { xdim = $2 }
/YDIMENSION/               { ydim = $2 }
/NROFLINES/                { nrows = $2 }
/NROFCELLSPERLINE/         { ncols = $2 }
/NROFBANDS/                { nbands = $2 }
/UNSIGNED8BITINTEGER/      { form = 0 }
/UNSIGNED16BITINTEGER/     { form = 1 }
/UNSIGNED32BITINTEGER/     { form = 2 }
/IEEE[48]BYTEREAL/         { print "ERROR - FLOAT NOT SUPPORTED BY GRASS" }
/ SIGNED8BITINTEGER/       { print "ERROR - SIGNED INTS NOT SUPPORTED BY GRASS" }
/HEADEROFFSET/             { print "ERROR - HEADEROFFSET NOT SUPPORTED BY GRASS - USE OUTPUT TO DATASET" }
/BYTEORDER/                { bytorder = $2 }
END {
            if (bytorder == LSBFIRST) print "ERROR - LSB byteorder NOT SUPPORTED BY GRASS";
            if (nbands > 1) print "ERROR - GRASS ONLY SUPPORTS ONE BAND - USE OUTPUT TO DATASET";
            if (substr(proj,1,7) ~ /\"*TM*/) print "proj:            1";
                 else print "proj:        UNKNOWN" ;
            if (substr(proj,1,7) ~ /\"*TM*/) zone = (substr(proj,8,2));
                 else zone = 1;
        printf("zone:            %s\n", zone);
        printf("north:           %f\n", regy + (offsety * ydim));
        printf("south:           %f\n", regy - ((nrows-offsety) * ydim));
        printf("east:            %f\n", regx + ((ncols-offsetx) * xdim));
        printf("west:            %f\n", regx - (offsetx * xdim));
        printf("cols:            %d\n", ncols);
        printf("rows:            %d\n", nrows);
        printf("e-w resol:       %f\n", xdim);
        printf("n-s resol:       %f\n", ydim);
        printf("format:          %d\n", form);
        printf("compressed:      0")
}
' >> $GRASS
-------------------------------------cut  here
#### END Import Shell Script ####

The script should be placed in your path ie./usr/local/bin
and executed with two file name arguments.
The first is the "absolute path" to the ERMapper image and the second is 
the name for the GRASS cell file
Please note I have not tried to use this script with LAT LONG datasets
so you might experience some problems here but it works well with UTM.

I have also included a script to take GRASS rasters to ERMapper which
you might also find useful.

#### Start Export Shell Script ####
-------------cut here
#!/bin/sh
#
# This shell script will create a .ERS file from a GRASS4.1 image header file 
# (cellhd directory) and link the cell image to the .ERS file location.
#

if [ "$#" -ne "2" ]
then
echo "Usage: export_grass_to_ermapper $LOCATION/cellhd/<file> $GISDBASE/$LOCATION_NAME/<file>"
exit 1
fi
GRASS=$LOCATION/cellhd/$1
ERS=$GISDBASE/$LOCATION_NAME/$2.ers
echo > $ERS "# ER Mapper $ERS file created from $GRASS file."
tr -d < $GRASS [:] | tr [a-z] [A-Z] | awk '
BEGIN {
        proj = 0;
        comp = 0;
        zone = 0;
        xdim = 0.0;
        ydim = 0.0;
        regx = 0.0;
        regy = 0.0;
        nbands = 1;
        ncols = 0;
        nrows = 0;
        nbits = 0;

}
/PROJ/                  { proj = $2 }
/ZONE/                  { zone = $2 }
/NORTH/                 { regy = $2 }
/WEST/                  { regx = $2 }
/COLS/                  { ncols = $2 }
/ROWS/                  { nrows = $2 }
/E-W/                   { xdim = $3 }
/N-S/                   { ydim = $3 }
/FORMAT/                { nbits = $2 }
/COMPRESSED/            { comp = $2 }
END {
print
printf("DatasetHeader Begin\n");
printf("\tVersion       = \"4.1\"\n");
printf("\tDataSetType   = ERStorage\n");
printf("\tDataType      = Raster\n");
printf("\tByteOrder     = MSBfirst\n");
printf("\tCoordinateSpace Begin\n");
printf("\tDatum = \"AGD66\"\n");
printf("\t%s%s%d%s\n","Projection       = ","\"TMAMG",zone,"\"");
printf("\tCoordinateType = EN\n");
printf("\tUnits = \"Meters\"\n");
printf("\tRotation      = 0:0:0.0\n");
printf("\tCoordinateSpace End\n");
printf("RasterInfo Begin\n");
        if(nbits=="0")
printf("\tCellType      = Unsigned8BitInteger\n");
        else
        if(nbits=="1")
printf("\tCellType      = Unsigned16BitInteger\n");
        else
        if(nbits=="2")
printf("\tCellType      = Unsigned32BitInteger\n");
printf("\tNullCellValue = 0\n");
printf("\tCellInfo Begin\n");
printf("\t\tXdimension = %f\n",xdim+0);
printf("\t\tYdimension = %f\n",ydim+0);
printf("\tCellInfo End\n");
printf("\tNrOfLines = %d\n",nrows+0);
printf("\tNrOfCellsPerLine = %d\n",ncols+0);
printf("\tRegistrationCoord Begin\n");
printf("\tEastings  = %d\n",regx+0);
printf("\tNorthings = %d\n",regy+0);
printf("\tRegistrationCoord End\n");
printf("\tNrOfBands = %d\n",nbands+0);
printf("\tRasterInfo End\n");
printf("DatasetHeader End\n");
}
' >> $ERS
ln -s $LOCATION/cell/$1 $GISDBASE/$LOCATION_NAME/$1
r.compress -u $1
-------------------------------------cut  here
#### END Export Shell Script ####

Note:ERMapper requires that the GRASS raster be uncompressed
     before it can read it.Therefore the above script runs
     r.compress on the cell file as its final task.This has
     some implications for management of disk storage space
     as some high resolution (small cell dimension) files can 
     become very large.

     Both scripts expect to be run from within GRASS as
     they rely on the GRASS environment variables.


        Good Luck
        Rod Paterson
        Aberfoyle Resources Ltd
        Phone: 61-3-8822226
        Fax:   61-3-8131086
        Email: rgp at aberfoyle.oz.au





More information about the grass-user mailing list