Issue with circular patterns showing up after orthographic projection

For a project I’m working on, I’m trying to project a GOES-16 image to an orthographic projection centred at the continental United States. In a later stage, I found these odd circular patterns appear in a result derived from an image I had reprojected to this orthographic projection. I have been able to trace this back to the projection, and can reproduce the patterns by computing the number of original GOES-16 pixels that are part of each pixel in the final orthographic reproduction.

Why is this? And if possible, how do I resolve it (besides choosing a different projection)?

I have included a (hopefully) working example below.

import numpy as np 
import cartopy.crs as ccrs
import pyproj
import pandas as pd
import matplotlib.pyplot as plt 


# Extent in orthographic projection coordinates
EXTENT = [-2373970.2054802035, 3638399.906890716, -2388159.0372323524, 1620933.528091551]  

# Orthographic projection
p_ortho = pyproj.Proj('+proj=ortho +ellps=GRS80 +lat_0=39.8283 +lon_0=-98.5795 +sweep=x +no_defs')

# GOES-16 nominal satellite height
h = 35786.0234375

# GOES-16 projection
proj_goes = pyproj.Proj(proj='geos', h=h*1000, lon_0=-75.2, sweep='x')

# GOES-16 scanning angles
x_goes = np.arange(-0.151844, 0.151844+5.5998564e-05, 5.5998564e-05)
y_goes = np.arange(0.151844, -0.151844-5.5998564e-05, -5.5998564e-05)


# Coordinates for lines of constant x and y scan angle 
xx_goes, yy_goes = np.meshgrid(x_goes[::100], y_goes)
lon_x, lat_x = proj_goes(xx_goes*h*1000, yy_goes*h*1000, inverse=True)

xx_goes, yy_goes = np.meshgrid(x_goes, y_goes[::100])
lon_y, lat_y = proj_goes(xx_goes*h*1000, yy_goes*h*1000, inverse=True)

# Create x-y grid to reproject
xx_goes, yy_goes = np.meshgrid(x_goes, y_goes)


# Project scanning angles to corresponding longitude latitude
lon_goes, lat_goes = proj_goes(xx_goes*h*1000, yy_goes*h*1000, inverse=True)

# Filter for CONUS domain (roughly)
lon = lon_goes[(lon_goes > -136)*(lon_goes < -45)*(lat_goes>10)*(lat_goes <55)]
lat = lat_goes[(lon_goes > -136)*(lon_goes < -45)*(lat_goes>10)*(lat_goes <55)]


# Now project longitudes, latitudes to orthographic projection
xx_p, yy_p = p_ortho(lon, lat)

# Use dataframe for calculation
df = pd.DataFrame({'x': xx_p, 'y' : yy_p})

# Min and max coordinates in orthographic projection
[x_min, x_max, y_min, y_max] = EXTENT

# Number of bins
x_bins = 300
y_bins = 200

# Calculate bin size
x_size = (x_max-x_min)/x_bins
y_size = (y_max-y_min)/y_bins

# Calculate bin number for xx_p, yyp
df['x_bin'] = np.floor(((df['x'] - x_min)/x_size)).astype(int)
df['y_bin'] = np.floor(((df['y'] - y_min)/y_size)).astype(int)

# Filter because I was lazy in defining the CONUS domain (?)
df = df[(df.x_bin >= 0) & (df.x_bin <=x_bins-1) & (df.y_bin >= 0) & (df.y_bin <=y_bins-1)]

# Reshape results into 2D array for display
res = np.zeros((y_bins, x_bins))
grouped = df.groupby(['x_bin', 'y_bin']).count().reset_index()
res[y_bins-1-grouped.y_bin.values, grouped.x_bin.values] = grouped.x.values

# Show results
fig = plt.figure(figsize=(10, 5), dpi=300)
proj = ccrs.Orthographic(central_latitude=39.8283, central_longitude=-98.5795)
ax = fig.add_axes([0, 0, 1, 1], projection=proj)
ax.set_extent(EXTENT, proj)
ax.set_aspect('auto')
cf = ax.imshow(res, origin="upper", extent=EXTENT)


# Lines of constant scan angle
for i in range(lat_x.shape[1]):
    
    ax.plot(lon_x[:,i], lat_x[:,i], linewidth=0.25, c='w', transform=ccrs.PlateCarree())
        
    ax.plot(lon_y[i,:], lat_y[i,:], linewidth=0.25, c='w', transform=ccrs.PlateCarree())
    
ax.coastlines()
fig.colorbar(cf, label="Number of GOES-16 ABI pixels per 100 orthographic projection pixels")