[MITgcm-support] Bathymetry file
Gus Correa
gus at ldeo.columbia.edu
Fri Mar 24 21:10:36 EDT 2017
Hi Jody, Saeid, Daniel
Maybe I am wrong, but the MITgcm, being written in Fortran,
should read and write in column-major order, not row-major.
Unless the developers twisted Fortran IO quite a bit,
which I don't think they did. [1]
Matlab, as Fortran, is column-major.
Hence, preparing binary files for the MITgcm in Matlab
should not require switching the multidimensional array
order from column-major to row-major, i.e. no transposition of 2D
arrays.
**
Your citation (from the MITgcm user guide, I think) says:
"...(x,y) ... ordered with the x coordinate varying fastest"
It goes like this (file or memory addresses increasing from the top):
[x,y] (x varying faster)
-----
(1,1)
(2,1)
(3,1)
...
(Nx,1) (Done with *column* 1 !)
(1,2)
(2,2)
(3,2)
...
(Nx,2) (Done with *column* 2!)
... etc
... etc
(1,Ny)
(2,Ny)
(3,Ny)
...
(Nx,Ny) (Done with *column* Ny)
So, the citation means that the data is in *column-major* format,
not row-major.
I.e., the first coordinate varies faster (the row index).
The first column comes first, then the second column, and so on.
**
Unfortunately, displaying 2D arrays
in Matlab is often the source of confusion [2].
Some commonly used Matlab display routines (e.g. contourf)
display 2D arrays, which Matlab seems consider as *matrices*,
with *columns* laid out *vertically*.
The first index (row index) increases
*from the bottom* to the top of the figure,
while the second index (column index)
increases from left to right.
That may be great for Linear-Algebra-oriented folks,
but not so much for those Geo- and Ocean-inclined.
Say, with Matlab's contourf(bathymetry), you get
a figure with this coordinate layout on the figure window:
(Nx,1)(Nx,2)..(Nx,Ny)
.....................
(2,1)(2,2)....(2,Ny)
(1,1)(1,2)....(1,Ny)
Hence, if your X represents longitude, and
Y latitude, you get a bathymetry "map" with longitudes varying
along the Y axis, and latitudes along X, which is rather weird.
North is to the left (not the top), South is to the right (not bottom).
The main ocean basins are laid out
horizontally, not vertically as one would expect.
Besides, this display is mirror-imaged, i.e.,
if you rotate the figure 90 degrees counterclockwise,
North and South go to the correct positions,
but the meridians are in decreasing order,
East is to the left and West to the right, i.e. flipped.
Please, see the contourf_bathymetry.pdf attached, made from the
verification/global_ocean.90x40x15 bathymetry.bin file,
with these Matlab commands:
>> [fid,msg]=fopen('bathymetry.bin','r','ieee-be');
>> bathy=fread(fid,40*90,'float32');
>> figure(1)
>> contourf(reshape(bathy,[90,40]))
>> colorbar
>> fclose(fid)
**
However, to display the 2D array as a *map*,
say, with X representing longitude and Y latitude,
South on the bottom, North on the top,
East to the right, West to the left,
which is what geo- or ocean-inclined people want to see,
one needs to transpose the 2D array for the figure to look right,
somehow giving it row-major *looks* (but it is *only the looks*!).
Say, with contourf(bathymetry')
(where the trailing single quote ' is the Matlab syntax for transpose)
you get a nice X=lon, Y=lat bathymetry map,
with this coordinate layout:
(1,Ny)(2,Ny)..(Nx,Ny)
.....................
(1,2)(2,2)....(Nx,2)
(1,1)(2,1)....(Nx,1)
Please, see the contourf_bathymetry_transposed.pdf attached,
made with these Matlab commands:
>> [fid,msg]=fopen('bathymetry.bin','r','ieee-be');
>> bathy=fread(fid,40*90,'float32');
>> figure(2)
>> contourf(reshape(bathy,[90,40])')
>> colorbar
>> fclose(fid)
**
Nevertheless, the need to transpose the array to display it as a lon-lat
map is a Matlab *display* issue only,
thanks to Matlab matrix-centered mindset.
It doesn't change the fact that in Matlab, in Fortran and
in MITgcm, the arrays are stored in memory and in binary files
in column-major order (the first coordinate varies faster).
My two cents,
Gus Correa
********************************************
[1] In the MDSIO code, files are opened as standard Fortran
direct access unformatted files, and read or written using standard
record-based Fortran READ and WRITE routines.
This means that array segments (i.e. file records,
often 2D-slices, but can be 3D-chunks too)
are read or written in Fortran column-major order
(i.e. first coordinate varies faster).
Please check the MDS_READ_FIELD and MDS_WRITE_FIELD routine codes:
http://mitgcm.org/public/r2_manual/latest/code_reference/vdb/code/1476.htm
and
http://mitgcm.org/public/r2_manual/latest/code_reference/vdb/code/1492.htm
*******
[2] Citing below from the Matlab contourf help page.
Note the X->column_index, Y->row_index
correspondence, and more specifically later,
the association of length(X) with size(Z,2),
and of length(Y) with size(Z,1).
Both imply that the *display* is transposed
(but not the underlying "Z" data array):
"contourf(Z) draws a filled contour plot of matrix Z, where Z is
interpreted as heights with respect to the x-y plane. Z must be at least
a 2-by-2 matrix that contains at least two different values. The x
values correspond to the column indices of Z and the y values correspond
to the row indices of Z. The contour levels are chosen automatically.
...
...
contourf(X,Y,Z), contourf(X,Y,Z,n), and contourf(X,Y,Z,v) draw filled
contour plots of Z using X and Y to determine the x and y values.
If X and Y are vectors, then length(X) must equal size(Z,2) and
length(Y) must equal size(Z,1). The vectors must be strictly increasing
or strictly decreasing and cannot contain any repeated values."
******
It can get even more confusing if to display the data
you use instead the Matlab routines such as imagesc,
because, as opposed to what contourf does,
imagesc lays out the array column indices
increasing from the top to the bottom of the figure.
(Argh, Matlab!)
That can be seen experimenting with contourf and imagesc
and the 2D array:
>> mm(:,1)=[1:1:20]
>> mm(:,2)=[1:1:20]
>> contourf(mm)
>> imagesc(mm)
On 03/23/2017 04:18 PM, Jody Klymak wrote:
> https://en.wikipedia.org/wiki/Row-_and_column-major_order
>
> So, matlab is column-major when it writes, numpy is row-major when it
> writes.
>
> The MITgcm reads in as row-major:
>
> "Line 46,
>
> bathyFile=’topog.box’
>
> This line specifies the name of the file from which the domain
> bathymetry is read. This file is a two-dimensional (x,y) map of depths.
> This file is assumed to contain 64-bit binary numbers giving the depth
> of the model at each grid cell, ordered with the x coordinate varying
> fastest. The points are ordered from low coordinate to high coordinate
> for both axes. The units and orientation of the depths in this file are
> the same as used in the MITgcm code. In this experiment, a depth of 0m
> indicates a solid wall and a depth of −5000m indicates open ocean. The
> matlab program input/gendata.m shows an example of how to generate a
> bathymetry file.”
>
>
>
>
>
>> On 23 Mar 2017, at 12:49 PM, Daniel Goldberg <dngoldberg at gmail.com
>> <mailto:dngoldberg at gmail.com>> wrote:
>>
>> Hi Saeid
>>
>> If you are using matlab, note that the ordering of the x and y
>> dimensions is opposite what is needed in your binary file. So if you
>> have a bathymetry array in matlab, you should do the transpose
>> operation before writing it to a binary file. (which is not the same
>> as rotating three times.)
>>
>> dan
>>
>> On Thu, Mar 23, 2017 at 7:03 PM, Menemenlis, Dimitris (329C)
>> <Dimitris.Menemenlis at jpl.nasa.gov
>> <mailto:Dimitris.Menemenlis at jpl.nasa.gov>> wrote:
>>
>> http://mitgcm.org/public/r2_manual/latest/online_documents/node102.html
>> <http://mitgcm.org/public/r2_manual/latest/online_documents/node102.html>
>>
>>
>>> On Mar 23, 2017, at 7:10 AM, smaeilpour <saeid_gal at yahoo.com
>>> <mailto:saeid_gal at yahoo.com>> wrote:
>>>
>>> Hi everyone,
>>> I have a question, to build a bathymetry file, is it necessary to
>>> rotate file 90 degrees and three times before i convert it to a
>>> binary file?
>>> any advice will be greatly appreciated.
>>>
>>> Regards,
>>> Saeid
>>
>>
>> _______________________________________________
>> MITgcm-support mailing list
>> MITgcm-support at mitgcm.org <mailto:MITgcm-support at mitgcm.org>
>> http://mitgcm.org/mailman/listinfo/mitgcm-support
>> <http://mitgcm.org/mailman/listinfo/mitgcm-support>
>>
>>
>>
>>
>> --
>>
>> Daniel Goldberg, PhD
>> Lecturer in Glaciology
>> School of Geosciences, University of Edinburgh
>> Geography Building, Drummond Street, Edinburgh EH8 9XP
>>
>>
>> em: D <mailto:dgoldber at mit.edu>an.Goldberg at ed.ac.uk
>> <mailto:an.Goldberg at ed.ac.uk>
>> web: http://ocean.mit.edu/~dgoldberg
>> _______________________________________________
>> MITgcm-support mailing list
>> MITgcm-support at mitgcm.org <mailto:MITgcm-support at mitgcm.org>
>> http://mitgcm.org/mailman/listinfo/mitgcm-support
>
> --
> Jody Klymak
> http://web.uvic.ca/~jklymak/
>
>
>
>
>
>
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: contourf_bathymetry.pdf
Type: application/pdf
Size: 282809 bytes
Desc: not available
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20170324/f6ab2046/attachment-0002.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: contourf_bathymetry_transposed.pdf
Type: application/pdf
Size: 276578 bytes
Desc: not available
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20170324/f6ab2046/attachment-0003.pdf>
More information about the MITgcm-support
mailing list