Issue with circular patterns showing up after orthographic projection #687

VincentMeijer · 2020-07-20T22:44:53Z

VincentMeijer
Jul 20, 2020

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.

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

I have included a (hopefully) working example below. Any help is much appreciated, and my apologies if this is not the correct place to ask questions like this.

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")