Two methods of displaying a starry sky. One generates a star map as a 64x64 png file using star position data and then uses that star map to display a sky shader. The other just displays random star positions along vertical strips.

bright_star_lst.csv

generate_star_skyamid.gd

@tool
extends EditorScript
# generates a panorama skyamid based on list of stars
# importing star data relies on CSV data importer add-on available here:
# https://github.com/timothyqiu/godot-csv-data-importer

var target_skyamid:String = "res://starry_sky/star_skyamid.png"

# bright star and near star csv files were derived from the spreadsheets here: 
# https://www.astrosci.scimuze.com/stellar_data.htm

const bright_star_data:String = "res://starry_sky/bright_star_lst.csv"

const near_star_data:String = "res://starry_sky/near_star_lst.csv"

# size of generated skyamid png
var target_size:int = 64

func pixel_to_vector(pixel:Vector2i) -> Vector3:
	var result:Vector3 = Vector3.ONE
	# convert pixel to range from -1.0 to 1.0 
	var scaled:Vector2 = ((Vector2(pixel) + Vector2.ONE * 0.5 ) / float(target_size) * 2.0 - Vector2.ONE)
	var rad:Vector2 = scaled * PI / 2.0
	if abs(scaled.x) > abs(scaled.y):
		if scaled.x < 0:
			scaled.y = -scaled.y
			rad.y = -rad.y
		var y:float = cos(rad.x)
		var dist:float = sin(rad.x)
		var side:float = scaled.y / abs(scaled.x) * PI / 4.0
		var x:float = cos(side) * dist
		var z:float = sin(side) * dist # * sign(scaled.y)
		#print(Vector3(x, y, z) , ", ", side, ", ", (Vector3(x, y, z).length()), " here")
		
		return Vector3(x, y, z)
	else:
		if scaled.y < 0:
			scaled.x = -scaled.x
			rad.x = -rad.x
		var y:float = cos(rad.y)
		var dist:float = sin(rad.y)
		var side:float = scaled.x / abs(scaled.y) * PI / 4.0
		var z:float = cos(side) * dist
		var x:float = sin(side) * dist # * sign(scaled.x)
		return Vector3(x, y, z)
	
	var angle_x:float = scaled.x 
	var angle_y:float = scaled.y
	
	result.x = (float(pixel.x) + 0.5) / float(target_size) * 2.0 - 1.0
	result.z = (float(pixel.y) + 0.5) / float(target_size) * 2.0 - 1.0
	result.y = sqrt(pow(result.x, 2.0) + pow(result.z, 2.0))
	return (result - Vector3(0.5, 0.5, 0.5)).normalized()

# include these just to test pixel_to_vector is working correctly
func test_inverse_pixel(eyedir:Vector3) -> Vector2i:
	var result:Vector2i = Vector2i( floor(test_inverse(eyedir) * float(target_size)))
	result.x = clamp(result.x, 0, target_size - 1)
	result.y = clamp(result.y, 0, target_size - 1)
	return result

func smooth_inverse_pixel(eyedir:Vector3) -> Vector2:
	var result:Vector2 = Vector2(test_inverse(eyedir) * float(target_size))
	return result

func test_inverse(eyedir:Vector3) -> Vector2:
	var ax:float = abs(eyedir.x);
	var xs:float = sign(eyedir.x);
	var az:float = abs(eyedir.z);
	var zs:float = sign(eyedir.z);
	var ay:float = abs(eyedir.y);
	var y:float = 1.0 - atan2(ay, Vector2(ax,az).length()) / 1.5708;
	var x:float = 0.0;
	if ax > az:
		x = xs * atan2(eyedir.z, ax) / 1.5708;
		return Vector2(0.5 + y * 0.5 * xs, 0.5 + y * x * xs);
	else:
		x = zs * atan2(eyedir.x, az) / 1.5708;
		return Vector2(0.5 + y * x * zs, 0.5 + y * 0.5 * zs);

# converts Right Ascension and Declanation to a direction vector 
func ra_dec_to_vector(ra:float, dec:float) -> Vector3:
	var result:Vector3
	# ra is in hours
	# dec is in degrees
	
	# convert to radians
	var ra_radians:float = ra / 24.0 * PI * 2.0
	var dec_radians:float = deg_to_rad(dec)
	
	var slice:Vector2 = Vector2(cos(dec_radians), sin(dec_radians))
	
	var side:Vector2 = Vector2(cos(-ra_radians), sin(-ra_radians)) * slice.x
	
	result.y = slice.y
	result.x = side.x
	result.z = side.y
	
	# TODO rotate for lat-long and time of day/year? 
	# convert to topocentric from geocentric?
	
	
	return result

func ra_to_float(ra:String) -> float:
	var split = ra.split_floats(" ")
	#print(split.size(), ", ", ra)
	return float(split[0]) + split[1]/60.0 * sign(split[0])

func dec_to_float(dec:String) -> float:
	var split = dec.split_floats(" ")
	#print(split.size(), ", ", dec)
	return float(split[0]) + split[1]/60.0 * sign(split[0])

func mag_to_alpha(mag:float) -> float:
	#var mag:float = float(mag_s)
	const max_mag:float = -1.44 # brightest mag we could have, will be equivalent to 1.0 alpha
	# for every 6 more than max_mag, we want 1/10th the alpha
	return 1.0 / pow(10.0, (mag - max_mag) / 6.0)

func _run() -> void:
	var star_positions = preload(bright_star_data)
	
	var edge_stars:Array[Vector4] = []
	var star_count:int = 0
	var skip_count:int = 0
	var star_image:Image = Image.create(target_size, target_size, false, Image.FORMAT_RGBA8)
	for i:int in range(0, star_positions.records.size()	):
		var ra:float = ra_to_float(star_positions.records[i].RA)
		var dec:float = dec_to_float(star_positions.records[i].DEC)
		var star_dir:Vector3 = ra_dec_to_vector(ra, dec)
		if star_dir.y < 0: # skip stars below horizon
			continue
		var pixel:Vector2i = test_inverse_pixel(star_dir)
		if star_image.get_pixel(pixel.x, pixel.y).g > 0:
			skip_count += 1
			continue
			
		var star_pixel_pos:Vector2 = smooth_inverse_pixel(star_dir)
		var offset:Vector2 = (star_pixel_pos - Vector2(pixel)) - Vector2(0.5, 0.5)
		if abs(offset.x) > 0.45 || abs(offset.y) > 0.45:
			print("found edge star")
			edge_stars.append(Vector4(star_dir.x, star_dir.y, star_dir.z, star_positions.records[i].V))
			
		star_count += 1
		# shift from -1 to 1 range to 0 to 1 range
		star_dir = star_dir * 0.5 + Vector3.ONE * 0.5
		var mag:float = mag_to_alpha(float(star_positions.records[i].V))
		var pixel_color:Color = Color(star_dir.x, star_dir.y, star_dir.z, mag)
		star_image.set_pixel(pixel.x, pixel.y, pixel_color)
	
	# do a pass to fill in adjacent pixels where the pixel is overlapping an edge
	for edge_star in edge_stars:
		var star_dir:Vector3 = Vector3(edge_star.x, edge_star.y, edge_star.z)
		var pixel:Vector2i = test_inverse_pixel(star_dir)
		var star_pixel_pos:Vector2 = smooth_inverse_pixel(star_dir)
		var offset:Vector2 = (star_pixel_pos - Vector2(pixel)) - Vector2(0.5, 0.5)
		var offset_pixel:Vector2i = pixel
		if abs(offset.x) > abs(offset.y):
			offset_pixel += Vector2i(sign(offset.x), 0)
		else:
			offset_pixel += Vector2i(0, sign(offset.y))
		if offset_pixel.x < 0 || offset_pixel.y < 0 || offset_pixel.x >= target_size || offset_pixel.y >= target_size:
			print("offset star out of bounds")
			continue
		if star_image.get_pixel(offset_pixel.x, offset_pixel.y).g > 0:
			print("skipping offset star")
			continue
		star_dir = star_dir * 0.5 + Vector3.ONE * 0.5
		var mag:float = mag_to_alpha(float(edge_star.w))
		var pixel_color:Color = Color(edge_star.x, edge_star.y, edge_star.z, mag)
		
		star_image.set_pixel(offset_pixel.x, offset_pixel.y, pixel_color)
	
	print(star_count, " ==> ", skip_count)
	
	var star_positions_02 = preload(near_star_data)
	
	# second list of stars, they had slightly different format
	
	for i:int in range(0, star_positions_02.records.size()	):
		var ra:float = (star_positions_02.records[i].RA)
		var dec:float = (star_positions_02.records[i].Dec)
		var star_dir:Vector3 = ra_dec_to_vector(ra, dec)
		if star_dir.y < 0: # skip stars below horizon
			continue
		var pixel:Vector2i = test_inverse_pixel(star_dir)
		if star_image.get_pixel(pixel.x, pixel.y).g > 0:
			skip_count += 1
			continue
			
		star_count += 1
		# shift from -1 to 1 range to 0 to 1 range
		star_dir = star_dir * 0.5 + Vector3.ONE * 0.5
		# adding a bonus column to let me cheat a few stars that I want to ensure are visible
		var mag:float = mag_to_alpha(float(star_positions_02.records[i].Vmag) + float(star_positions_02.records[i].Bonus))
		var pixel_color:Color = Color(star_dir.x, star_dir.y, star_dir.z, mag)
		star_image.set_pixel(pixel.x, pixel.y, pixel_color)
	
	print(star_count, ", ", skip_count)
	
	star_image.save_png(target_skyamid)

skyamid_starry_sky.gdshader

shader_type sky;
// displays stars on skybox
//

uniform vec4 background_gradient_horizon:source_color = vec4(0.0, 0.0, 0.15, 1.0);
uniform vec4 background_gradient_top:source_color = vec4(0.0, 0.0, 0.0, 1.0);

// range of sizes for star
// basically in millionths?
uniform vec2 star_size_range = vec2(2, 10);

// texture that defines position and strength of every star
uniform sampler2D skyTexture:filter_nearest;

vec2 get_skyamid_uv(vec3 eyedir) {
	float ax = abs(eyedir.x);
	float xs = sign(eyedir.x);
	float az = abs(eyedir.z);
	float zs = sign(eyedir.z);
	float ay = abs(eyedir.y);
	float y = 1.0 - atan(ay, length(vec2(ax,az))) / 1.5708;
	float x = 0.0;
	if (ax > az) {
		x = xs * atan(eyedir.z, ax) / 1.5708;
		return vec2(0.5 + y * 0.5 * xs, 0.5 + y * x * xs);
	}
	else {
		x = zs * atan(eyedir.x, az) / 1.5708;
		return vec2(0.5 + y * x * zs, 0.5 + y * 0.5 * zs);
	}
}

void sky() {
	float adjusted_y = clamp(pow(EYEDIR.y, 1.0), 0.0 ,1.0);
	vec3 background = mix(background_gradient_horizon.rgb, background_gradient_top.rgb, adjusted_y);

	// pick pixel from skyamid
	vec2 sky_uv = get_skyamid_uv(EYEDIR);
	//sky_uv = round(sky_uv * 64.00);
	vec4 sample = texture(skyTexture, sky_uv);
	// use that for direction of star
	float star_strength = mix(1.0 - star_size_range.x * 0.000001, 1.0 - star_size_range.y * 0.000001, sample.a);
	
	vec3 star_dir = normalize(sample.xyz * 2.0 - vec3(1.0)); // assumes direction is directly
	float aligned = smoothstep(star_strength, 1.00000001, dot(EYEDIR, star_dir));

	COLOR = background + aligned * sqrt(sample.a);
}

starry_sky.gdshader

shader_type sky;
// displays stars as sky shader
// one star per vertical strip but with enough strips it passes first inspection

uniform vec4 background_gradient_horizon:source_color = vec4(0.0, 0.0, 0.15, 1.0);
uniform vec4 background_gradient_top:source_color = vec4(0.0, 0.0, 0.0, 1.0);
// number of vertical strips with one star in each
uniform float vertical_strips = 300.0;
// rangel of angle
uniform float vertical_scale:hint_range(0.0, 1.57, 0.01) = 1.57;
// power on vertical to adjust distribution
uniform float vertical_pow = 2.0;

// range of sizes for star
// basically in millionths?
uniform vec2 star_size_range = vec2(2, 10);

float random_height(float seed) {
	return fract(tan(seed * 120.3) * 137.0 + tan(seed * 12.0) * 31.0);
}

void sky() {
	float adjusted_y = clamp(pow(EYEDIR.y, 1.0), 0.0 ,1.0);
	vec3 background = mix(background_gradient_horizon.rgb, background_gradient_top.rgb, adjusted_y);
	// find vertical strip
	float rounded_angle = round(atan(EYEDIR.x, EYEDIR.z) * vertical_strips / 2.0 / PI) / vertical_strips * 2.0 * PI;
	// angle height of star in this strip
	float star_angle = pow(random_height(rounded_angle), vertical_pow) * vertical_scale;
	float star_strength = mix(1.0 - star_size_range.x * 0.000001, 1.0 - star_size_range.y * 0.000001, random_height(rounded_angle + 3.1));
	vec3 star_dir = normalize(vec3(sin(rounded_angle) * cos(star_angle), sin(star_angle), cos(rounded_angle) * cos(star_angle)));
	float aligned = smoothstep(star_strength, 1.0, dot(EYEDIR, star_dir));
	COLOR = background + aligned;
}