implement lat_lon_to_l_c equivalent to goesr.l1b.PugFile

himawari/ahi2cmi.py

@@ -41,6 +41,12 @@ SCENE_ID_CONVERSION = {
     "FLD": "Full Disk",
 }
 
+# FUTURE: if available, fetch base class
+# try:
+#     from goesr.l1b import PugFile
+# except ImportError:
+#     PugFile = object
+
 class HimawariAHIasCMIP(object):
     """
     yield path, value, dimensions, attributes frames used by goesr.cmi_changer
@@ -420,6 +426,7 @@ class HimawariAHIasCMIP(object):
             yield self.p('%s_image_bounds' % name), vrange, self.d(DEFAULT_YX_BOUNDS_NAME), a
 
         # y and x variables and their bounds and centerpoints
+        # negate scaling of y to conform to GOES / PROJ4 convention
         for each in yx_frames('y', y, mb.my, mb.by, negate_scaling=True):
             yield each
 
@@ -490,6 +497,34 @@ class HimawariAHIasCMIP(object):
         bwl = np.array([bwl_um], dtype=np.float32)
         yield self.p('band_id'), band, self.d(DEFAULT_BAND_DIM_NAME), a
 
+    def lat_lon_to_l_c(self, lat, lon, round=True):
+        """
+        convert latitude-longitude all the way back to line-column by reversing projection and y/x scaling
+        :param lat: latitude value or array in degrees
+        :param lon: longitude value or array in degrees
+        :return: line, column array indexed 0..height-1, 0..width-1
+        """
+        from pyproj import Proj
+        proj = Proj(projparams=dict(self.proj4_params))
+        # convert to nadir-meters using proj4
+        xm, ym = proj(lon, lat, inverse=False)
+        # divide out height to get back to projection view angles
+        h = self.perspective_point_height
+        y, x = ym / h, xm / h
+        _, _, scale = self._hs.coords(lines=1, columns=1, unscaled=True)
+        # extract scale factor and add offset to back out to line-from-0, column-from-0
+        # also note Y is inverse between AHI and GOES/PROJ4
+        ysf, yao = -scale.my, -scale.by
+        xsf, xao = scale.mx, scale.bx
+        # reverse scale factor and add offset to get line and column
+        # note that AHI coordinate system counts from 1! so we subtract 1.0 to get consistency with GOES
+        line, column = (y - yao) / ysf - 1.0, (x - xao) / xsf - 1.0
+        # round to nearest integer
+        if round:
+            return np.round(line).astype(np.int32), np.round(column).astype(np.int32)
+        else:  # return floating point values
+            return line, column
+
     def walk(self, extra_data_attrs={}):
         """
         iterate PVDA frames equivalent to a PUG nc_walk