Problem subsetting this GRIB2 dataset

[GoogleCodeExporter]

geosubetting doesn't work for this GRIB2 dataset in NCTOOLBOX (But it works in 
the older NJTBX).  Here's the details, reported by Dave Thompson 
(dave.thompson3 at gmail.com), and confirmed by me:

I'm attempting to subset the data found in any of the four data files found 
here: 
http://weather.noaa.gov/pub/SL.us008001/ST.expr/DF.gr2/DC.ndgd/GT.slosh/AR.conus
/

I used to do this with njToolBox as:
>>nc = mDataset('grib2.mdlsurgegrid.00con');
>>lon = [-82.18 -80.18];
>>lat = [26.52 30.52];
>>[data,grd] =  nj_subsetGrid(nc,'Extra_Tropical_Storm_Surge',[lon,lat]);
>> grd

grd = 

     lat: [92x48 double]
     lon: [92x48 double]
    time: [97x1 double]
       z: []

>> size(data)

ans =

         97.00         92.00         48.00

With nctoolbox:
>> nc = ncgeodataset('grib2.mdlsurgegrid.00con');
log4j:WARN No appenders could be found for logger (ucar.nc2.NetcdfFile).
log4j:WARN Please initialize the log4j system properly.
>> s.lon = [-82.18 -80.18];
>> s.lat = [26.52 30.52];
>> gvar = nc.geovariable('Extra_Tropical_Storm_Surge');
>> sub = gvar.geosubset(s)
Error using ncdataset/findvariable (line 436)
Java exception occurred:
java.lang.NullPointerException

at ucar.unidata.util.EscapeStrings.backslashEscape(EscapeStrings.java:554)

at ucar.nc2.NetcdfFile.escapeName(NetcdfFile.java:278)

Error in ncdataset/attributes (line 257)
                v = obj.findvariable(variable);

Error in ncgeodataset/geovariable (line 147)
                att = obj.attributes(variableName);

Error in ncgeodataset/subsref (line 367)
                            B = builtin('subsref',obj,g);

Error in ncgeovariable/getlonvar (line 347)
            lv = src.dataset.geovariable(tn);

Error in ncgeovariable/getlondata (line 411)
            v = src.getlonvar;

Error in ncgeovariable/geoij (line 758)
                    g.lon = obj.getlondata([1, 1, 1], obj.size, [1, 1, 1]);

Error in ncgeovariable/geosubset (line 566)
            [indstart_r indend_r indstart_c indend_c] = obj.geoij(struct);

Error in ncgeovariable/subsref (line 969)
                            sref = builtin('subsref',obj,s);

I can see that the grib file's coordinate variables are 'x' & 'y' which are km 
in Lambert Conformal projection rather than 'lon' and 'lat'. Is this the 
problem for nctoobox?

Thank you,
Dave

Original issue reported on code.google.com by rsignell on 5 Nov 2012 at 2:56

[GoogleCodeExporter]

Looks like a problem in java with escaping the variable name. I'm going to 
assign to Brian since he implemented the variable escaping code.

Original comment by crosb...@gmail.com on 5 Nov 2012 at 3:35

[GoogleCodeExporter]

Starting investigation ... here's the CDL for grib2.mdlsurgegrid.00con:

netcdf /Users/brian/Downloads/grib2.mdlsurgegrid.00con {
 dimensions:
   time = 97;
   y = 689;
   x = 1073;
 variables:
   float Extra_Tropical_Storm_Surge(time=97, y=689, x=1073);
     :units = "m";
     :long_name = "Extra_Tropical_Storm_Surge @ surface";
     :missing_value = NaNf; // float
     :grid_mapping = "Lambert_Conformal";
     :GRIB_param_discipline = "Oceanographic_products";
     :GRIB_param_category = "Surface_Properties";
     :GRIB_param_name = "Extra_tropical_storm_surge";
     :GRIB_generating_process_type = "Forecast";
     :GRIB_param_id = 2, 10, 3, 193; // int
     :GRIB_product_definition_template = 0; // int
     :GRIB_product_definition_template_desc = "Analysis/forecast at horizontal level/layer at a point in time";
     :GRIB_level_type = 1; // int
     :GRIB_level_type_name = "surface";
     :GRIB_VectorComponentFlag = "easterlyNortherlyRelative";
   char Lambert_Conformal;
     :grid_mapping_name = "lambert_conformal_conic";
     :standard_parallel = 25.0; // double
     :longitude_of_central_meridian = 265.0; // double
     :latitude_of_projection_origin = 25.0; // double
     :earth_shape = "Earth spherical with radius specified by producer in m";
     :earth_radius = 6371200.0; // double
     :GRIB_param_Dx = 5079.406; // double
     :GRIB_param_Dy = 5079.406; // double
     :GRIB_param_GDSkey = 801839731; // int
     :GRIB_param_La1 = 20.192; // double
     :GRIB_param_LaD = 25.0; // double
     :GRIB_param_Latin1 = 25.0; // double
     :GRIB_param_Latin2 = 25.0; // double
     :GRIB_param_Lo1 = 238.446; // double
     :GRIB_param_LoV = 265.0; // double
     :GRIB_param_NpProj = "true";
     :GRIB_param_Nx = 1073; // int
     :GRIB_param_Ny = 689; // int
     :GRIB_param_ProjFlag = 0; // int
     :GRIB_param_Quasi = "false";
     :GRIB_param_ResCompFlag = 0; // int
     :GRIB_param_SpLat = -90.0; // double
     :GRIB_param_SpLon = 0.0; // double
     :GRIB_param_VectorComponentFlag = "easterlyNortherlyRelative";
     :GRIB_param_Winds = "True";
     :GRIB_param_grid_name = "Lambert_Conformal";
     :GRIB_param_grid_radius_spherical_earth = 6371200.0; // double
     :GRIB_param_grid_shape = "Earth spherical with radius specified by producer in m";
     :GRIB_param_grid_shape_code = 1; // int
     :GRIB_param_grid_type = 30; // int
     :GRIB_param_grid_units = "m";
     :GRIB_param_scanning_mode = 64; // int
     :_CoordinateTransformType = "Projection";
     :_CoordinateAxisTypes = "GeoX GeoY";
   int time(time=97);
     :long_name = "forecast time";
     :units = "hour since 2012-11-05T00:00:00Z";
     :_CoordinateAxisType = "Time";
   double y(y=689);
     :units = "km";
     :long_name = "y coordinate of projection";
     :standard_name = "projection_y_coordinate";
     :grid_spacing = "5.07940576171875 km";
     :_CoordinateAxisType = "GeoY";
   double x(x=1073);
     :units = "km";
     :long_name = "x coordinate of projection";
     :standard_name = "projection_x_coordinate";
     :grid_spacing = "5.07940576171875 km";
     :_CoordinateAxisType = "GeoX";

 :Conventions = "CF-1.4";
 :Originating_center = "US National Weather Service - NCEP(WMC) (7)";
 :Originating_subcenter = "NWS Meteorological Development Laboratory(MDL) (14)";
 :Generating_Model = "Probabilistic Storm Surge";
 :Product_Status = "Operational test products";
 :Product_Type = "Forecast products";
 :title = "US National Weather Service - NCEP(WMC) Probabilistic Storm Surge Forecast products";
 :institution = "Center US National Weather Service - NCEP(WMC) (7) Subcenter NWS Meteorological Development Laboratory(MDL) (14)";
 :source = "Type: Forecast products Status: Operational test products";
 :history = "Direct read of GRIB-2 into NetCDF-Java 4 API";
 :CF:feature_type = "GRID";
 :file_format = "GRIB-2";
 :location = "/Users/brian/Downloads/grib2.mdlsurgegrid.00con";
 :_CoordinateModelRunDate = "2012-11-05T00:00:00Z";
}

Original comment by bschlin...@gmail.com on 5 Nov 2012 at 5:38

[GoogleCodeExporter]

Found the error:

In ncgeovariable.m,  there is a method named 'getlonname'. The contents of the 
function are:

function ln = getlonname(src)
    % NCGEOVARIABLE.getlonname()
    for i = 1:length(src.axes)
        tempname = src.axes{i};
        javaaxisvar  =   src.dataset.netcdf.findVariable(tempname);
        type{i} = char(javaaxisvar.getAxisType());
    end
    match = strcmp('Lon', type);
    ln = src.axes(match);
end

The key point to note is that it's looking for 'Lon' to match against. However 
the axis types returned for this file are:

type = 

    'Time'    'GeoY'    'GeoX'

I'll modify the gelonname and getlatname functions to accept either Lat or GeoY 
and Lon or GeoX. Give me a few minutes to make and commit the changers.

Original comment by bschlin...@gmail.com on 5 Nov 2012 at 5:54

[GoogleCodeExporter]

OK, I'm bouncing this back to Alex. Alex, If I make the following changes:

diff -r 9b3debfeb51a cdm/ncgeovariable.m
--- a/cdm/ncgeovariable.m   Wed Jul 25 11:27:53 2012 -0400
+++ b/cdm/ncgeovariable.m   Mon Nov 05 09:59:37 2012 -0800
@@ -392,7 +392,7 @@
                 javaaxisvar  =   src.dataset.netcdf.findVariable(tempname);
                 type{i} = char(javaaxisvar.getAxisType());
             end
-            match = strcmp('Lon', type);
+            match = strcmp('Lon', type) | strcmp('GeoX', type);
             ln = src.axes(match);
         end

@@ -403,7 +403,7 @@
                 javaaxisvar  =   src.dataset.netcdf.findVariable(tempname);
                 type{i} = char(javaaxisvar.getAxisType());
             end
-            match = strcmp('Lat', type);
+            match = strcmp('Lat', type) | strcmp('GeoY', type);
             ln = src.axes(match);
         end

Then a different exception is thrown by geoij. The stack trace is:

Error using ncgeovariable/geoij (line 809)
Longitude contains values that follow both -180/180 and 0/360+ conventions; can 
not subset.

Error in ncgeovariable/geosubset (line 599)
            [indstart_r indend_r indstart_c indend_c] = obj.geoij(struct);

Error in ncgeovariable/subsref (line 1002)
                            sref = builtin('subsref',obj,s);

So, the issue is that the data file uses 'km' as units instead of degrees. 
geoij will need to be modified to deal with north and east geovalues that are 
not in degrees and fall outside the +/-180 or 0-360 value ranges.

Original comment by bschlin...@gmail.com on 5 Nov 2012 at 6:02

[GoogleCodeExporter]

Commited my changes listed above as a 'partial fix'. 

Original comment by bschlin...@gmail.com on 5 Nov 2012 at 6:08

[GoogleCodeExporter]

Interesting, Our code should be doing the conversion from projected coordinates 
to lat/lon based on identification as GeoY or GeoX, and thats not happening. 
Ill have to take a look at at code with the newest changes.

Original comment by crosb...@gmail.com on 5 Nov 2012 at 6:10

[GoogleCodeExporter]

I must have created this bug myself when i refactored

Original comment by crosb...@gmail.com on 5 Nov 2012 at 6:15

[GoogleCodeExporter]

I have committed changes to fix the problem. Please verify that these work. You 
can get the latest code from either this google code hg repository or you can 
download a package of the latest unstable dev code at 
https://github.com/acrosby/nctoolbox_recent/archive/master.zip

Original comment by crosb...@gmail.com on 5 Nov 2012 at 8:57

• Changed state: Fixed

[GoogleCodeExporter]

Works! Thank you!

Original comment by dave.tho...@gmail.com on 6 Nov 2012 at 12:34

[GoogleCodeExporter]

Hey Alex, Rich, Is there any reason we shouldn't roll and release a new version 
of nctoolbox with this fix?

Original comment by bschlin...@gmail.com on 6 Nov 2012 at 5:01

[GoogleCodeExporter]

No reason that I can think of

Original comment by crosb...@gmail.com on 6 Nov 2012 at 6:52

[GoogleCodeExporter]

OK, new release has been posted to downloads.

Original comment by bschlin...@gmail.com on 6 Nov 2012 at 7:04