[GRASSLIST:892] Re: Problems importing a GMT grd file (SOLUTION)
David Finlayson
dfinlays at u.washington.edu
Sun Aug 3 13:06:48 EDT 2003
This is a wrap-up to a problem I had importing GMT grd files into GRASS
using r.in.bin on Cygwin/XP.
Apparently, GMT 3.4.X type=1 grids are platform dependent. The header
line in type=1 grids created on Windows (Cygwin or precompiled binaries)
contains 4 extra bytes that are not found in a type=1 grid created on
Linux. (for more information search the GMT-HELP archives for
[GMT-HELP:11231] Problems with grdreformat on Cygwin and Windows )
r.in.bin assumes that the header in type=1 GMT grids conforms to the
Linux convention, hence, any type=1 grids created on the Windows
platform will fail to import properly.
The solution provided here is a small program that can strip the extra 4
bytes out of the Windows type=1 grids and convert them to the Linux
convention (or vis-à-vis). Thanks to Marta E. Ghidella for this program.
How to proceed:
1) Create the GMT grid in Cygwin or DOS version of GMT as usual (Note
the default grid type is type=0)
2) Use GMT's grdreformat to convert the grid to type=1:
> grdreformat ingrid=0 outgrid.win=1
3) Run convert6432.exe or convert6432.py as follows:
> ./convert6432.exe outgrid.win outgrid.lin -4
or alternatively
> python convert6432.py outgrid.win outgrid.lin -4
4) Start grass and use r.in.bin to import the linux-style GMT grid:
GRASS:~> r.in.bin -hf in=outgrid.lin out=grid
That, at lest, works for me on Windows XP with Cygwin
NOTE: To compile the C version of convert6432.exe in Cygwin, download
convert6432.c to any cygwin directory and type:
$ gcc -mno-cygwin convert6432.exe -o convert6432.c
There is a fixed memory size in the C version that is adjusted by making
the variable NMAX=4000000 larger and recompiling. This is unnecessary
in the Python version so long as you don't run out of swap space.
--
David Finlayson
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
-------------- next part --------------
/*
convert6432.c
Program to modify the header in grd type 1 GMT files
If these files are generated in Linux, the size of the header is 892;
if they are created in Cygwin or by the Microsoft C binaries, the size
is 896.
This program modifies the header, thus converting files from one
implementation to the other.
>From Linux to Cygwin, it inserts 4 bytes in the header
>From Cygwin to Linux, it removes 4 bytes in the header
The original file is not modified; a new converted file is created.
convert6432 <input grd=1 file> <output grd=1 file> ishift
ishift = 4 to add 4 bytes (linux to cygwin)
ishift = -4 to remove 4 bytes (cygwin to linux)
Marta, August 3, 2003
mghidella at dna.gov.ar
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <conio.h>
#include <io.h>
main (argc,argv)
int argc;
char **argv;
{
char *infile, *outfile, *bufch, chaux[4];
FILE *fpin,*fpout;
/* Note that NMAX can be larger; it has to be greater than or equal to the
grid size in bytes */
int num,ie,imem,NMAX=4000000, iread,iw,i;
int nx,ny,node_offset,ishift,npad=0,*pnx;
float *grid;
if (argc != 4) {
fprintf(stderr,"Needed argumets:\n");
fprintf(stderr,"1) Input grd file (grd format =1)\n");
fprintf(stderr,"2) Output grd file (grd format =1)\n");
fprintf(stderr,"3) Number of bytes to shift in the header:\n");
fprintf(stderr," 4 (adds 4 bytes, linux to cygwin conversion)\n");
fprintf(stderr," -4 (removes 4 bytes, cygwin to linux conversion)\n");
exit(-1);
}
infile=argv[1]; outfile=argv[2];
ishift=atoi(argv[3]);
printf("infile: %s; outfile: %s; shift= %d\n",infile,outfile,ishift);
if (ishift != 4 && ishift != -4) {
fprintf(stderr,"third argument has to be either 4 or -4;\n exits ...\n");
exit(-1);
}
if(fpin=fopen(infile,"rb"))
fprintf(stderr, "opened file: %s \n", infile);
else {fprintf(stderr, "error opening: %s \n", infile);exit(-1);}
if(fpout=fopen(outfile,"wb"))
fprintf(stderr, "opened file: %s \n", outfile);
else {fprintf(stderr, "error opening: %s \n", outfile);exit(-1);}
/* Reads the first 4 bytes in the file to check that it is not a netcdf file */
fread(chaux,1,4,fpin);
if (chaux[0]==67 && chaux[1]==68 && chaux[2]==70 && chaux[3]==1) {
fprintf(stderr, "seems to be a netcdf file: chaux= %s \nexits...\n",chaux);
exit(-1);
}
/* gets nx from auxiliary char */
pnx = (int *) chaux;
nx=*pnx;
/* Reads ny and node_offset */
fread(&ny,4,1,fpin);
fread(&node_offset,4,1,fpin);
/* writes ny, ny and node_offset in the output file */
fwrite(&nx,4,1,fpout);
fwrite(&ny,4,1,fpout);
fwrite(&node_offset,4,1,fpout);
printf("nx= %d; ny= %d, node_offset=%d\n",nx,ny,node_offset);
/* num is the numbr of grid elements */
num=nx*ny;
/* only reads in this case */
if(ishift == -4) {fread(&npad,4,1,fpin);
}
/* only writes in this case */
if(ishift == 4) {fwrite(&npad,4,1,fpout);
}
printf("num= %d,num; npad= %d\n",num,npad);
/* calculates number of bytes needed for storage */
imem=sizeof(float)*num;
/* checks memory */
if(imem > NMAX) {
fprintf(stderr,"needed memory greater than expected ...");
exit(-1);
}
else
fprintf(stderr,"imem= %d\n",imem);
/* allocates memory for the remainder of the header */
bufch=(char *) malloc(880);
/* allocates grid memory */
grid=(float *) malloc(imem);
/* reads and writes the rest of the header */
iread=fread(bufch,1,880,fpin);
iw=fwrite(bufch,1,880,fpout);
printf("iread= %d iw= %d\n",iread,iw);
/* reads the whole grid */
iread=fread(grid,4,num,fpin);
/* checks for errors */
if( iread != num ) {
fprintf(stderr,"conversion failed; exits ...\n");
exit(-1);
}
/* writes the entire grid */
iw=fwrite(grid,4,num,fpout);
printf("iread= %d iw= %d\n",iread,iw);
free(grid);
free(bufch);
fclose(fpin);
fclose(fpout);
}
-------------- next part --------------
# convert6432.py
#
# Program to modify the header in grd type 1 GMT files
# If these files are generated in Linux, the size of the header is 892;
# if they are created in Cygwin or by the Microsoft C binaries, the size
# is 896.
# This program modifies the header, thus converting files from one
# implementation to the other.
# From Linux to Cygwin, it inserts 4 bytes in the header
# From Cygwin to Linux, it removes 4 bytes in the header
# The original file is not modified; a new converted file is created.
#
# convert6432 <input grd=1 file> <output grd=1 file> ishift
#
# ishift = 4 to add 4 bytes (linux to cygwin)
# ishift = -4 to remove 4 bytes (cygwin to linux)
#
# Marta, August 3, 2003
#
# mghidella at dna.gov.ar
#
# Translated into Python by David Finlayson (dfinlays at u.washington.edu)
import sys
import struct
# Get Command line Arguments
if len(sys.argv) != 4:
print """
USAGE: python convert6432.py <input grd=1 file> <output grd=1 file> ishift
Needed arguments:
1) Input grd file (grd format =1)
2) Output grd file (grd format = 1)
3) Number of bytes to shift in the header:
4 (adds 4 bytes, linux to cygwin conversion)
-4 (removes 4 bytes, cygwin to linux conversion)
"""
sys.exit(-1)
infile = sys.argv[1]
outfile = sys.argv[2]
ishift = int(sys.argv[3])
print "infile: %s; outfile: %s; shift= %d" % (infile, outfile, ishift)
if (ishift != 4 and ishift != -4):
print "third argument has to be either 4 or -4"
sys.exit(-1)
try:
fin = file(infile, 'rb')
except:
print "error opening: %s" % (infile)
sys.exit(-1)
try:
fout = file(outfile, 'wb')
except:
print "error opening: %s" % (outfile)
sys.exit(-1)
# Reads the first 4 bytes in the file to check that it is not a netCDF file
aux = fin.read(4)
auxstr = []
auxstr.append(struct.unpack("c", aux[0])[0])
auxstr.append(struct.unpack("c", aux[1])[0])
auxstr.append(struct.unpack("c", aux[2])[0])
auxstr.append(str(struct.unpack("b", aux[3])[0]))
auxstr = "".join(auxstr)
if auxstr == "CDF1":
print "seems to be a netCDF file: auxstr= %s" % (auxstr)
sys.exit(-1)
# gets nx from auxilary char
nx = struct.unpack("i", aux)[0]
# Reads ny and node_offset
ny = fin.read(4)
ny = struct.unpack("i", ny)[0]
node_offset = fin.read(4)
node_offset = struct.unpack("i", node_offset)[0]
# writes nx, ny and node_offset in the output file
fout.write(struct.pack("i", nx))
fout.write(struct.pack("i", ny))
fout.write(struct.pack("i", node_offset))
print "nx= %d; ny= %d, node_offset=%d" % (nx, ny, node_offset)
# only reads in this case
if ishift == -4:
npad = fin.read(4)
# only writes in this case
if ishift == 4:
fout.write(struct.pack("i", 0))
# reads and writes the rest of the file
fout.write(fin.read())
fin.close()
fout.close()
More information about the grass-user
mailing list