Merge pull request #160 from dopplershift/gini-fixes

GINI fixes

metpy/io/cdm.py

@@ -379,11 +379,17 @@ def cf_to_proj(var):
         kwargs['lon_0'] = var.longitude_of_central_meridian
         kwargs['lat_1'] = var.standard_parallel
         kwargs['lat_2'] = var.standard_parallel
-    if var.grid_mapping_name == 'polar_stereographic':
+    elif var.grid_mapping_name == 'polar_stereographic':
         kwargs['proj'] = 'stere'
         kwargs['lon_0'] = var.longitude_of_projection_origin
         kwargs['k_0'] = 1.0    # scale factor at natural origin
         kwargs['x_0'] = False  # Easting
         kwargs['y_0'] = False  # Northing
+    elif var.grid_mapping_name == 'mercator':
+        kwargs['proj'] = 'merc'
+        kwargs['lon_0'] = var.longitude_of_projection_origin
+        kwargs['lat_ts'] = var.standard_parallel
+        kwargs['x_0'] = False  # Easting
+        kwargs['y_0'] = False  # Northing
 
     return pyproj.Proj(**kwargs)

metpy/io/gini.py

@@ -257,48 +257,28 @@ def to_dataset(self):
             proj_var.longitude_of_central_meridian = self.proj_info.lov
             proj_var.latitude_of_projection_origin = self.prod_desc2.lat_in
             proj_var.earth_radius = 6371200.0
-            proj = cf_to_proj(proj_var)
+            _add_projection_coords(ds, self.prod_desc, proj_var, self.proj_info.dx,
+                                   self.proj_info.dy)
         elif self.prod_desc.projection == GiniProjection.polar_stereographic:
             proj_var = ds.createVariable('Polar_Stereographic', np.int32)
             proj_var.grid_mapping_name = 'polar_stereographic'
-            proj_var.standard_parallel = self.prod_desc2.lat_in
             proj_var.longitude_of_projection_origin = self.proj_info.lov
-            proj_var.latitude_of_projection_origin = self.prod_desc2.lat_in
+            proj_var.latitude_of_projection_origin = -90 if self.proj_info.proj_center else 90
+            proj_var.earth_radius = 6371200.0
+            _add_projection_coords(ds, self.prod_desc, proj_var, self.proj_info.dx,
+                                   self.proj_info.dy)
+        elif self.prod_desc.projection == GiniProjection.mercator:
+            proj_var = ds.createVariable('Mercator', np.int32)
+            proj_var.grid_mapping_name = 'mercator'
+            proj_var.longitude_of_projection_origin = self.prod_desc.lo1
+            proj_var.latitude_of_projection_origin = self.prod_desc.la1
+            proj_var.standard_parallel = self.prod_desc2.lat_in
             proj_var.earth_radius = 6371200.0
-            proj = cf_to_proj(proj_var)
+            _add_projection_coords(ds, self.prod_desc, proj_var, self.prod_desc2.resolution,
+                                   self.prod_desc2.resolution)
         else:
             raise NotImplementedError('Need to add more projections to dataset!')
 
-        # Get projected location of lower left point
-        x0, y0 = proj(self.prod_desc.lo1, self.prod_desc.la1)
-
-        # Coordinate variable for x
-        ds.createDimension('x', self.prod_desc.nx)
-        x_var = ds.createVariable('x', np.float64, dimensions=('x',))
-        x_var.units = 'm'
-        x_var.long_name = 'x coordinate of projection'
-        x_var.standard_name = 'projection_x_coordinate'
-        x_var[:] = x0 + np.arange(self.prod_desc.nx) * (1000. * self.proj_info.dx)
-
-        # Now y
-        ds.createDimension('y', self.prod_desc.ny)
-        y_var = ds.createVariable('y', np.float64, dimensions=('y',))
-        y_var.units = 'm'
-        y_var.long_name = 'y coordinate of projection'
-        y_var.standard_name = 'projection_y_coordinate'
-        y_var[:] = y0 + np.arange(self.prod_desc.ny) * (1000. * self.proj_info.dy)
-
-        # Get the two-D lon,lat grid as well
-        x, y = np.meshgrid(x_var[:], y_var[:])
-        lon, lat = proj(x, y, inverse=True)
-        lon_var = ds.createVariable('lon', np.float64, dimensions=('y', 'x'), wrap_array=lon)
-        lon_var.long_name = 'longitude'
-        lon_var.units = 'degrees_east'
-
-        lat_var = ds.createVariable('lat', np.float64, dimensions=('y', 'x'), wrap_array=lat)
-        lat_var.long_name = 'latitude'
-        lat_var.units = 'degrees_north'
-
         # Now the data
         name = self.prod_desc.channel
         if '(' in name:
@@ -332,3 +312,38 @@ def __str__(self):
                  'Lower Left Corner (Lon, Lat): ({0.lo1}, {0.la1})',
                  'Resolution: {1.resolution}km']
         return '\n\t'.join(parts).format(self.prod_desc, self.prod_desc2)
+
+
+def _add_projection_coords(ds, prod_desc, proj_var, dx, dy):
+    'Helper function for adding coordinate variables (projection and lon/lat) to a dataset'
+    proj = cf_to_proj(proj_var)
+
+    # Get projected location of lower left point
+    x0, y0 = proj(prod_desc.lo1, prod_desc.la1)
+
+    # Coordinate variable for x
+    ds.createDimension('x', prod_desc.nx)
+    x_var = ds.createVariable('x', np.float64, dimensions=('x',))
+    x_var.units = 'm'
+    x_var.long_name = 'x coordinate of projection'
+    x_var.standard_name = 'projection_x_coordinate'
+    x_var[:] = x0 + np.arange(prod_desc.nx) * (1000. * dx)
+
+    # Now y
+    ds.createDimension('y', prod_desc.ny)
+    y_var = ds.createVariable('y', np.float64, dimensions=('y',))
+    y_var.units = 'm'
+    y_var.long_name = 'y coordinate of projection'
+    y_var.standard_name = 'projection_y_coordinate'
+    y_var[:] = y0 + np.arange(prod_desc.ny) * (1000. * dy)
+
+    # Get the two-D lon,lat grid as well
+    x, y = np.meshgrid(x_var[:], y_var[:])
+    lon, lat = proj(x, y, inverse=True)
+    lon_var = ds.createVariable('lon', np.float64, dimensions=('y', 'x'), wrap_array=lon)
+    lon_var.long_name = 'longitude'
+    lon_var.units = 'degrees_east'
+
+    lat_var = ds.createVariable('lat', np.float64, dimensions=('y', 'x'), wrap_array=lat)
+    lat_var.long_name = 'latitude'
+    lat_var.units = 'degrees_north'

metpy/io/tests/test_gini.py

@@ -23,7 +23,7 @@ def test_gini_basic():
     assert pdb.channel == 'WV (6.5/6.7 micron)'
     assert pdb.num_records == 1280
     assert pdb.record_len == 1100
-    assert pdb.datetime, datetime(2015, 12, 8, 22, 0, 19 == 0)
+    assert pdb.datetime == datetime(2015, 12, 8, 22, 0, 19)
     assert pdb.projection == GiniProjection.lambert_conformal
     assert pdb.nx == 1100
     assert pdb.ny == 1280
@@ -58,7 +58,7 @@ def test_gini_ak_regional():
     assert pdb.channel == 'IR (3.9 micron)'
     assert pdb.num_records == 408
     assert pdb.record_len == 576
-    assert pdb.datetime, datetime(2016, 4, 8, 14, 45, 20 == 0)
+    assert pdb.datetime == datetime(2016, 4, 8, 14, 45, 20)
     assert pdb.projection == GiniProjection.polar_stereographic
     assert pdb.nx == 576
     assert pdb.ny == 408
@@ -84,6 +84,43 @@ def test_gini_ak_regional():
     assert f.data.shape, (pdb.num_records == pdb.record_len)
 
 
+def test_gini_mercator():
+    'Test reading of GINI file with Mercator projection (from HI)'
+    f = GiniFile(get_test_data('HI-REGIONAL_4km_3.9_20160616_1715.gini'))
+    pdb = f.prod_desc
+    assert pdb.source == 1
+    assert pdb.creating_entity == 'GOES-15'
+    assert pdb.sector_id == 'Hawaii Regional'
+    assert pdb.channel == 'IR (3.9 micron)'
+    assert pdb.num_records == 520
+    assert pdb.record_len == 560
+    assert pdb.datetime == datetime(2016, 6, 16, 17, 15, 18)
+
+    assert pdb.projection == GiniProjection.mercator
+    assert pdb.nx == 560
+    assert pdb.ny == 520
+    assert_almost_equal(pdb.la1, 9.343, 4)
+    assert_almost_equal(pdb.lo1, -167.315, 4)
+
+    proj = f.proj_info
+    assert proj.resolution == 0
+    assert_almost_equal(proj.la2, 28.0922, 4)
+    assert_almost_equal(proj.lo2, -145.878, 4)
+    assert proj.di == 0
+    assert proj.dj == 0
+
+    pdb2 = f.prod_desc2
+    assert pdb2.scanning_mode == [False, False, False]
+    assert_almost_equal(pdb2.lat_in, 20.0, 4)
+    assert pdb2.resolution == 4
+    assert pdb2.compression == 0
+    assert pdb2.version == 1
+    assert pdb2.pdb_size == 512
+    assert pdb2.nav_cal == 0
+
+    assert f.data.shape, (pdb.num_records == pdb.record_len)
+
+
 def test_gini_bad_size():
     'Test reading a GINI file that reports a bad header size'
     f = GiniFile(get_test_data('NHEM-MULTICOMP_1km_IR_20151208_2100.gini'))
@@ -115,6 +152,31 @@ def test_gini_ak_regional_dataset():
     assert ds.variables['IR'].grid_mapping in ds.variables
     assert_almost_equal(ds.variables['lon'][0, 0], f.prod_desc.lo1, 4)
     assert_almost_equal(ds.variables['lat'][0, 0], f.prod_desc.la1, 4)
+    assert_almost_equal(ds.variables['Polar_Stereographic'].longitude_of_projection_origin,
+                        210.0)
+    assert_almost_equal(ds.variables['Polar_Stereographic'].latitude_of_projection_origin,
+                        90.0)
+
+
+def test_gini_mercator_dataset():
+    'Test the dataset interface for a GINI file with Mercator projection'
+    f = GiniFile(get_test_data('HI-REGIONAL_4km_3.9_20160616_1715.gini'))
+    ds = f.to_dataset()
+    assert 'x' in ds.variables
+    assert 'y' in ds.variables
+    assert 'IR' in ds.variables
+    assert hasattr(ds.variables['IR'], 'grid_mapping')
+    assert ds.variables['IR'].grid_mapping in ds.variables
+    lat = ds.variables['lat']
+    lon = ds.variables['lon']
+    assert_almost_equal(lon[0, 0], f.prod_desc.lo1, 4)
+    assert_almost_equal(lat[0, 0], f.prod_desc.la1, 4)
+    # 2nd corner lat/lon are at the "upper right" corner of the pixel, so need to add one
+    # more grid point
+    assert_almost_equal(lon[-1, -1] + (lon[-1, -1] - lon[-1, -2]), f.proj_info.lo2, 4)
+    assert_almost_equal(lat[-1, -1] + (lat[-1, -1] - lat[-2, -1]), f.proj_info.la2, 4)
+    assert_almost_equal(ds.variables['Mercator'].longitude_of_projection_origin, -167.315)
+    assert_almost_equal(ds.variables['Mercator'].latitude_of_projection_origin, 9.343)
 
 
 def test_gini_str():