Coastline

[henryb]

Hi there,
I have some NOAA Wave Watch data at 10 min and 30 min resolution and am trying to get a nice coastline with MCIDAS-V
More specifically I am plotting Significant wave height (HTSGW)

With FERRET I went thou the following procedure to get a nice coastine
1) Create a sea mask from from a 5min topographic nc file (etopo05.nc)
2) Fill missing values in HTSGW with average nearest neighbour values in lon and lat.
( this is a custom ferret function -fill_XY() )
3) Reinterpolate HTSGW on the 5min grid, apply sea mask , and then call contour/fill function

Should I go through the same procedure with McIDAS (if so how ?) or can I do something clever with the "World Coastlines" shape file ?
I have noticed a few functions that take a shapefile map as an argument:
subsetFromMap()
subsetRangeFromMap()

Attached is a image of Florida HTSGW and the Coastline and Legend

Many thanks

[Rick]

Hi Henry,

If I understand your question correctly, your land sea mask looks jagged and you'd like to clean up the graphic - correct? This might be just a matter of placing the mask on top of the wave data. If the jagged appearance is from the wave data not reaching the land, we might have to find a way of filling in the data. If you could, please add your netCDF files them as attachment to the forum. If they are too large, please post them to our anonymous ftp location.

ftp ftp.ssec.wisc.edu
login: anonymous
password: e-mail address
cd pub
cd incoming
bin
put {files}

[henryb]

HI Rick,
thank you very much for your assistance.
The jagged appearance is from wave data - I haven't actually applied the mask
I would just like to fill the area from wave data to coastline shape with sensible values.
Ideally I would like a method that works at different zoom levels - The NOAA file I have
are global - 30m, atlantic -10 min, atlantic - 4m


I have posted the at_10m netcdf file :
just the first time frame: the variables lon,lat,HTSGW

Thanks
... Henry

[Rick]

Hi Henry,

One of programmers gave me the piece of code (at the bottom) and it looks like it might meet your needs. Copy the code below and paste it into your jython library. From the main display ...
Tools>Formulas>Jython Library

Now create a formula to use the new code.
Tools>Formulas>Create Formula

You can add anything for the name. For the formula, enter the following:
smooth_coast(wave_height)

From the Data Explorer load your data from the General File Chooser.

From the Field Selector choose formulas and then the name of the formula created. For the example I attached, I chose color shaded plan view.

Code: Select all

def smooth_coast(a):
  import Numeric as np;
  b = a[0].getFloats();
  c = np.reshape(b[0], [331,301]);
  for i in xrange(1,330):
    for j in xrange(1,300):
      if c[i][j]!=c[i][j]:
         cnt = 0;
         sum = 0;
         for ii in xrange(i-1,i+2):
           for jj in xrange(j-1,j+2):
              if c[ii][jj] == c[ii][jj]:
                 cnt = cnt + 1;
                 sum = sum + c[ii][jj];
         if cnt > 3:
           c[i][j] = sum / cnt;
  b[0] = np.reshape(c, [331*301]);
  a[0].setSamples(b);
  return a;

The programmer also noted that you might want to fiddle with the "cnt > 3" line, you might get some better results.

[henryb]

HI Rick,
thank you very much for your help.
I've given the formula a quick go - it kind of works
Im a little confused by the line:
if c[i][j]!=c[i][j]:


... Henry

[henryb]

HI rick,
OK so the line
c[i]!=c[i]
check for NaN
cool

... Henry

[tomw]

HI Rick,
thank you very much for your help.
I've given the formula a quick go - it kind of works
Im a little confused by the line:
if c[i][j]!=c[i][j]:

Henry --

I know it looks weird, but that is a test for a missing value (Not-a-Number, or NaN). Jython has no explicit way to do that, so I used this form. By definition, if a value is an IEEE NaN then it will not be equal to iteself.

You might want to change the value in the "cnt > 3" line. Also, please note that I just used the fixed dimensions from the sample data; if you change the data source, you will need to change those....if that's the case, you can pick up the dimensions using:

Code: Select all

dims = b.getDomainSizes()

which will return an int[] of the lengths of the two dimensions.

tom

[henryb]

Hi Rick & tom,
In the end I managed to get a nice coastline by
1) replacing the missing values in HTSGW (sig wave height) by 0.0
2) Overlaying a topology mask which is 0 -for z <0 and 1 for z >=0.
In addition the mask is transparent for z <0

The problem I now have is that the 1 minute topographic mask over the region of interest is eating loads of memory.

1) cut down the size of the mask by breaking it up into small regions
2) Would it be possible to create the mask from the default world coastlines ?
I've tried to fill the map without any real success


... Henry

[Rick]

Hi Henry,

What is the file format of you topo map? Some data formats, you can specify a stride which samples the data.

Rick

[henryb]

Hi Rick,
The file format is netcdf -
so the 1min map is in dims - 21601x10801
5min map is in dims - 4320x2161


Even when I spatially sub-set the 1min map for the region of interest e.g
the gulf of mexico Itswallowing up 4G of memory