Computing desk
< January 21	<< Dec \| January \| Feb >>	January 23 >

Welcome to the Wikipedia Computing Reference Desk Archives
The page you are currently viewing is a transcluded archive page. While you can leave answers for any questions shown below, please ask new questions on one of the current reference desk pages.

January 22

Strange result from code to make an arbitrary azimuthal equidistant projection from an equirectangular projection

Hi Python coders,

I wish to render Danny Quah's Valeriepieris circle in azimuthal equidistant projection centred on a given latitude and longitude, by transforming the equirectangular projection File:Earthmap1000x500.jpg, and found formulas for the transformation from (x, y) to (lat, long) at http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html. However my Python 2 code below which uses PyPNG gives me the image in the thumbnail.

Can someone please point out where it went wrong (and a possible remedy)?

Thanks,
cmɢʟee⎆τaʟκ 01:51, 22 January 2022 (UTC)[reply]

#!/usr/bin/env python

path_in       = 'mya/Earthmap1000x500.png'
path_out      = 'Valeriepieris_circle_azimuthal_equidistant.png'
lat_centre    = 21.7
long_centre   = 99.383333
out_size      = 1024
out_size_half = out_size // 2

class Png:
 def __init__(self, path_in):
  (self.width, self.height, self.pixels, self.metadata) = png.Reader(path_in).read_flat()
  self.planes = self.metadata['planes']
 def __str__(self): return str((self.width, self.height, len(self.pixels), self.metadata))
 def write(self, path_out):
  png.Writer(width=self.width, height=self.height,
             bitdepth=self.metadata['bitdepth'], interlace=self.metadata['interlace'],
             planes=self.metadata['planes'], greyscale=self.metadata['greyscale'],
             alpha=self.metadata['alpha']).write_array(open(path_out, 'wb'), self.pixels)

import re, math, png

## Formulas from http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html
def azimuthal_equidistant_to_equirectangular(x, y, lat1, long1):
 c        = math.radians(math.hypot(x, y))
 if c == 0 or (abs(lat1) == 90 and y == 0): return (0, 0)
 sin_c    = math.sin(c)
 cos_c    = math.cos(c)
 lat1_rad = math.radians(lat1)
 sin_lat1 = math.sin(lat1_rad)
 cos_lat1 = math.cos(lat1_rad)
 to_asin  = cos_c * sin_lat1 + y * sin_c * cos_lat1 / c
 if abs(to_asin) > 1: return (0, 0)
 lat  =  math.degrees(math.asin(to_asin))
 long = (math.degrees(math.atan2(-x, y) if lat1 ==  90 else
                      math.atan2( x, y) if lat1 == -90 else
                      math.atan2(x * sin_c, c * cos_lat1 * cos_c - y * sin_lat1 * sin_c)) +
         long1 + 540) % 360 - 180
 # print(x, y, lat, long)
 return (lat, long)

png_in  = Png(path_in)
print(png_in)
print(png_in.pixels[:20])
png_out = Png(path_in) ## copy most of original's metadata
png_out.width  = png_out.height = out_size
png_out.pixels = [0] * (png_out.width * png_out.height)
print(png_out)
for  out_y in range(out_size):
 for out_x in range(out_size):
  (lat, long) = azimuthal_equidistant_to_equirectangular(
                 out_x - out_size_half, out_size_half - out_y, lat_centre, long_centre)
  in_y = int(png_in.height * ( 90 - lat ) / 180.0)
  in_x = int(png_in.width  * (180 + long) / 360.0)
  in_offset  = (in_y  * png_in.width + in_x ) * png_in .planes
  out_offset = (out_y * out_size     + out_x) * png_out.planes
  png_out.pixels[out_offset:out_offset + png_out.planes] = (
   png_in.pixels[in_offset :in_offset  + png_in.planes ])
 # print(y, in_offset, out_offset)
png_out.write(path_out)

I would suggest asking this question on a more targeted Q&A site such as Stack Overflow, where you may find more viewers with the necessary skills to help. That said, my first impression is that the issue may be one of units. Your code is working in a mix of pixel coordinates, degrees, and radians. I suspect you are missing a conversion somewhere. In particular, your out_x/out_y variables appear to be pixel coordinates, but the first two parameters of azimuthal_equidistant_to_equirectangular seem to be expecting an angle (degrees?), presumably limited to some expected range. -- Tom N talk/contrib 17:02, 23 January 2022 (UTC)[reply]

A superficial examination of the formulas at MathWorld suggests that the values for the first two parameters of azimuthal_equidistant_to_equirectangular should be given with the Earth radius as unit, so the value of math.hypot(x, y) should not exceed 1.0. --Lambiam 23:37, 23 January 2022 (UTC)[reply]

Many thanks! You're right. I've now fixed it as above. Cheers, cmɢʟee⎆τaʟκ 01:45, 23 February 2022 (UTC) [reply]

Resolved

Get Firefox to open files in new tab in last active window

Using Firefox on Mac – if I have the browser set as my default to open, say, jpegs from Finder, how do I get it to stop opening them in new windows, and to instead open them in a new tab in the last active window? I've set browser.link.open_newwindow.override.external to 3, but it doesn't work. --Lazar Taxon (talk) 17:27, 22 January 2022 (UTC)[reply]

Linux command-line program for JPG analysis

I am looking for a command-line program for Linux that would take a JPG file as an input parameter, and instead of displaying it, would print information about it on stdout. The only real requirement is that it should print the image dimensions, i.e. width and height. Is there one? JIP | Talk 20:09, 22 January 2022 (UTC)[reply]

The file command will print basic information about a jpg file. If you need more detailed information you might look at ImageMagick. CodeTalker (talk) 20:18, 22 January 2022 (UTC)[reply]

exiftool will give you all the metadata. You can filter using the tool, or else pipe the output into either grep or awk which is less efficient but more flexible.

 $ exiftool O_come_O_come.jpg | grep Image
 Image Width                     : 6392
 Image Height                    : 4240
 Image Size                      : 6392x4240

for example. Martin of Sheffield (talk) 21:20, 22 January 2022 (UTC)[reply]

If you have the ImageMagick package installed you can also use identify:

    $ identify O_come_O_come.jpg
    O_come_O_come.jpg JPEG 6392x4240 6392x4240+0+0 8-bit sRGB 424627B 0.000u 0:00.000

This also works on several other formats, such as PNG and GIF. --Lambiam 23:15, 23 January 2022 (UTC)[reply]

Thanks Lambian. The actual output is O_come_O_come.jpg JPEG 6392x4240 6392x4240+0+0 8-bit Gray 256c 4490570B 0.000u 0:00.000 Martin of Sheffield (talk) 23:21, 23 January 2022 (UTC)[reply]