diff --git a/metpy/io/cdm.py b/metpy/io/cdm.py index b940e3bbbb0..9a1d76a5b10 100644 --- a/metpy/io/cdm.py +++ b/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) diff --git a/metpy/io/gini.py b/metpy/io/gini.py index 658c78ded49..7cb67fd6af7 100644 --- a/metpy/io/gini.py +++ b/metpy/io/gini.py @@ -257,47 +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.longitude_of_projection_origin = self.proj_info.lov proj_var.latitude_of_projection_origin = -90 if self.proj_info.proj_center else 90 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.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 + _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: @@ -331,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' diff --git a/metpy/io/tests/test_gini.py b/metpy/io/tests/test_gini.py index 6e79785f637..a1eac31684b 100644 --- a/metpy/io/tests/test_gini.py +++ b/metpy/io/tests/test_gini.py @@ -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')) @@ -121,6 +158,27 @@ def test_gini_ak_regional_dataset(): 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(): 'Test the str representation of GiniFile' f = GiniFile(get_test_data('WEST-CONUS_4km_WV_20151208_2200.gini')) diff --git a/testdata/HI-REGIONAL_4km_3.9_20160616_1715.gini b/testdata/HI-REGIONAL_4km_3.9_20160616_1715.gini new file mode 100644 index 00000000000..94c8a2b7736 Binary files /dev/null and b/testdata/HI-REGIONAL_4km_3.9_20160616_1715.gini differ