Creates a haversine version of a raster elevation file for use in Vue Infinite.

create-haversine-raster.py

class HaversineRaster(object):
    """Creates Haversine elevation raster for input into Vue
    
    """
    def __init__(self, project_dir, project_name, raster):
        self.project_dir = project_dir
        self.project_name = project_name
        self.in_raster = raster
        self.props = v_utils.get_raster_props(raster)
        self.hav_raster = self.create_haversine_raster()
        self.hav_props = v_utils.get_raster_props(self.hav_raster)
    
    def create_haversine_raster(self):
        """Calculate a Haversine distance raster"""
        ras_to_pts = os.path.join(self.project_dir, self.project_name + ".gdb", 
                                  self.project_name + "_RasterToPoints")
        hav_dist_raster = os.path.join(self.project_dir, 
                                       self.project_name + ".gdb", 
                                       self.project_name + "_HavDist")
        
        msg = "Converting raster to points, please be patient, "\
              "as this will take a while..."
        arcpy.AddMessage(msg)
        logging.info(msg)
        arcpy.RasterToPoint_conversion(self.in_raster, ras_to_pts, "VALUE")
        
        new_fields = [["Hav_Distance", "Float", "0"]]
        v_utils.add_fields(ras_to_pts, new_fields)
        
        # Get POINT_X and POINT_Y fields with lat/lons
        arcpy.AddXY_management(ras_to_pts)  
        
        lat = self.props["y_center"]
        lon = self.props["x_center"]
        expression = "haversine(!POINT_X!, !POINT_Y!)"
        # Codeblock to calculate haversine distance in meters into point FC
        codeblock = """def haversine(lon2, lat2):
                            lon1 = {0}
                            lat1 = {1}
                            lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
                            dlon = lon2 - lon1
                            dlat = lat2 - lat1
                            a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
                            c = 2 * math.asin(math.sqrt(a))
                            m = {2} * c
                            return m""".format(lon, lat, v_utils.EARTH_RADIUS)    
        arcpy.CalculateField_management(ras_to_pts, "Hav_Distance", 
                                        expression, "PYTHON", codeblock)
        # Create raster with haversine distance from raster center
        msg = "Creating new haversine distance raster: " + self.project_name + "_HavDist..."
        arcpy.AddMessage(msg)
        logging.info(msg)
        arcpy.PointToRaster_conversion(ras_to_pts, "Hav_Distance", hav_dist_raster,
                                       "MOST_FREQUENT", "Hav_Distance", self.in_raster)
        arcpy.env.extent = self.in_raster
        self.out_haversine_raster_gcs = os.path.join(self.project_dir, self.project_name + ".gdb",
                                           self.project_name + "_HavCurvedGCS")
        # Raster calculation to get Haversine elevation raster
        msg = "Creating haversine raster for Vue..."
        arcpy.AddMessage(msg)
        logging.info(msg)
        self.out_curved_hav_raster_gcs = arcpy.Raster(self.in_raster) - ((arcpy.Raster(hav_dist_raster) * arcpy.Raster(hav_dist_raster)) / (2 * v_utils.EARTH_RADIUS))
        self.out_curved_hav_raster_gcs.save(self.out_haversine_raster_gcs)
        msg = "GCS haversine raster computed: " + self.out_haversine_raster_gcs
        arcpy.AddMessage(msg)
        logging.info(msg)
        # Project curved elevation raster
        # Get center of GCS original raster
        gcs_raster_props = v_utils.get_raster_props(self.in_raster)
        gcs_x_center = gcs_raster_props["x_center"]
        gcs_y_center = gcs_raster_props["y_center"]
        # Create custom spatial reference, putting raster center x/y as long/lat of center
        custom_sr = 'PROJCS["VueGnomonic",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",\
                     SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],\
                     UNIT["Degree",0.0174532925199433]],PROJECTION["Gnomonic"],\
                     PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],\
                     PARAMETER["Longitude_Of_Center",{0}],PARAMETER["Latitude_Of_Center",{1}],\
                     UNIT["Meter",1.0]]'.format(gcs_x_center, gcs_y_center)
        self.out_haversine_raster_pcs = os.path.join(self.project_dir, 
                                                     self.project_name + ".gdb",
                                                     self.project_name + "_HavCurvedPCS")
        msg = "Projecting haversine raster: " + self.out_haversine_raster_pcs + "..."
        arcpy.AddMessage(msg)
        logging.info(msg)
        # No geog transformation required since the raster we use to calculate
        #   haversine on is already in WGS84 datum
        arcpy.ProjectRaster_management(self.out_curved_hav_raster_gcs, 
                                       self.out_haversine_raster_pcs,
                                       custom_sr, "BILINEAR")
        return self.out_haversine_raster_pcs