Good evening everyone, many thanks for previous answers to my question in this forum, I have a FORTRAN program that is used to project geographic grid data to rectangular grid form. I tried to use the FORTRAN program but it didn’t work, then I read the help page of the FORTRAN program which says my geographic grid must be in standard USGS grid form before it can work. The FORTRAN code is part of a Potential field software suite that works only on Windows XP.
My question is how can I convert my geographic grid data ( grid file format = netCDF format) to the standard USGS grid format?
Thanks
The GDAL package has a program called “gdal_translate” that can be used to convert raster files to different formats. I am not sure what your Fortran program means by “standard USGS grid” format, but one of the formats that is supported by “gdal_translate” is called “USGSDEM”. This is an old ASCII file format that the USGS used to use for their DEM products. The GDAL description is “USGS Optional ASCII DEM (and CDED)”. You could use a command like:
gdal_translate -of USGSDEM DEM_gmt.grd DEM_USGS.txt
to translate your data to the USGSDEM format.
You can do that directly in GMT with the grdconvert module. But if the intention is to reproductive the grid, why not use grdproject directly?
reproject
Damn correctors
I did not know about the grdconvert command. I looked at the manpage, and I am not sure if it uses GDAL for its conversion or not.
I understood that @John needed the output in a specific old USGS format, so the question was about the file format, not the projection.
All modules, not only grdconvert, that write grids can save them through GDAL. See this.
I mentioned grdproject because @John wanted to project a grid and was referring to a fortran program to do that. So if USGSDEM is the driver name that gdal_translate would take, then the grdproject would for example be something like this.
grdproject input.grd -J... -Gresult.grd=gd:USGSDEM
@EJFielding and @Joaquim many many thanks, am so grateful. @EJFielding I actually didn’t get the way you wrote the script. But the issue is I have a grid file, just the normal grid files like free-air gravity and Etopo1 grid and I tried to project them at a particular grid spacing but the result returned variable grid spacing rather than integer grid spacing then I got the Potential field software, I saw the Fortran code that can project from geographic coordinate to rectangular coordinate but I tried to use the normal grid, but the code gave an error and requested I transform my data from original grid form to USGS standard grid before I can use the code. Basically this code will give me the gird spacing in integer rather than variable. And I have tried many things here but didn’t work out. So @Joaquim was on point. MY MAIN AIM ON THIS GRID IS JUST TO ACHIEVE AN INTEGER GRID SPACING RATHER THAN VARIABLE. I wish to upload the FORTRAN code but the attached only allows image files.
thanks
I expect that GMT6 and GMT5 work the same way as GMT4 for the grid spacing. If you want an exact integer grid spacing you can add “=” to the end of the “-D” spacing flag. Then it adjusts the region to make the grid spacing exact. For example, grdproject -D20k= input.grd -Goutput.grd should make the grid spacing exactly 20 km (of course you also need to specify the region and projection).
@EJFielding, many thanks to you, am so grateful. Yes, I got it this time with your script. I am so grateful, I just used -D20= with other flags, and it worked well. But how come if a person try to use Mercator projection such that -Jm/1:1 and -Jm/1:1 don’t give same X_max and Y_max for each projection used.
Or which is better to use for Mercator projection purpose? note I have tried this but want more clarity.
many thanks to all, grateful