API tools faq
paste
Naeddyr

Equirectangular rotation, G'mic / GIMP / Krita [updated Jan 2026: 'endif' > 'fi', fix float divby0]

Naeddyr
Aug 23rd, 2020 (edited)
983
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 60.49 KB | None | 0 0
raw download clone embed print report
  1. #@gui ____<b>Map Projection</b>
  2. #-----------------------------
  3. ###############################
  4. ## EQUIRECTANGULAR ROTATION
  5. #@gmic rotate_equirectangular_map : roll angle, pitch angle, yaw angle, revert bool
  6. #@gmic : Take a map or equirectangular panorama
  7. #@gmic : and rotate it around three axises
  8.  
  9. #@gui Rotate Equirectangular Map : rotate_equirectangular_map
  10. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order in degrees.")
  11. #@gui : Roll = float(0,-360,360)
  12. #@gui : Pitch ("north and south") = float(0,-360,360)
  13. #@gui : Yaw ("east and west") = float(0,-360,360)
  14. #@gui : Revert rotations = bool(false)
  15. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  16.  
  17.  
  18. rotate_equirectangular_map :
  19.  
  20. # Default to 0 degrees
  21. -skip ${1=0},${2=0},${3=0},${4=0}
  22.  
  23. # Black magic spell to make the preview look better???
  24. _fF_zoomaccurate=1
  25. _fF_zoomfactor=1
  26.  
  27.  
  28. # Create new image, the ""-block is the mathematical 'function' that
  29. # determines the value of the output pixel wherever we are iterating
  30.  
  31. -input 100%,100%,1,100%,"
  32.  
  33. # Command parameters from the command line or the sliders.
  34. # Variables from outside the ""-block have to be outside the ""s.
  35. yaw="{$3}";
  36. pitch="{$2}";
  37. roll="{$1}";
  38. revert="{$4}";
  39.  
  40. # Converting degrees into radians for the matrices later
  41. # Those minus signs are just a direction correction.
  42. alpha=-roll * pi / 180;
  43. beta=-pitch * pi / 180;
  44. gamma=-yaw * pi / 180;
  45.  
  46. # Equirectangular coordinates into spherical radians coordinates
  47. # phi <- height
  48. # theta <- width
  49. # radius rho=1, we are always working with a unit sphere / globe with radius 1
  50. # The extra pi in "otheta=PI..." is to offset the axis towards the center of the image
  51. otheta=pi + (x/w)*2*pi;
  52. ophi=(y/h)*pi;
  53. rho=1;
  54.  
  55. # Convert to Cartesian x, y, z, we can leave out rho=1
  56. # cx=rho*sin(ophi)*cos(otheta);
  57. # cy=rho*sin(ophi)*sin(otheta);
  58. # cz=rho*cos(ophi);
  59.  
  60. cx=sin(ophi)*cos(otheta);
  61. cy=sin(ophi)*sin(otheta);
  62. cz=cos(ophi);
  63.  
  64. # Do the transformation in Cartesian using rotation matrix maths magic
  65. if( revert==0,
  66.  
  67. rotz=cz*( cos(alpha)*cos(beta) )
  68. + cy*( cos(alpha)*sin(beta)*sin(gamma) - sin(alpha)*cos(gamma) )
  69. + cx*( cos(alpha)*sin(beta)*cos(gamma) + sin(alpha)*sin(gamma) );
  70.  
  71. roty=cz*( sin(alpha)*cos(beta) )
  72. + cy*( sin(alpha)*sin(beta)*sin(gamma) + cos(alpha)*cos(gamma) )
  73. + cx*( sin(alpha)*sin(beta)*cos(gamma) - cos(alpha)*sin(gamma) );
  74.  
  75. rotx=cz*(-sin(beta))
  76. + cy*(cos(beta)*sin(gamma))
  77. + cx*(cos(beta)*cos(gamma)); ,
  78.  
  79. # !!!!!!!!!!!!!!!! Reversion matrix??
  80. # small letter x = cos(x)
  81. # capital letter X = sin(x)
  82. # {{ a*b, -A*b, B }
  83. # { B*G*a+A*g, -A*B*G+a*g, -G*b }
  84. # { A*G-B*a*g, G*a+A*B*g, b*g }}
  85.  
  86. ## Testing reversion
  87.  
  88. rotz=cz*( cos(-alpha)*cos(-beta) )
  89. + cy*( -sin(-alpha)*cos(-beta) )
  90. + cx*( sin(-beta) );
  91.  
  92. roty=cz*( sin(-beta)*sin(-gamma)*cos(-alpha) + sin(-alpha)*cos(-gamma) )
  93. + cy*( -sin(-alpha)*sin(-beta)*sin(-gamma) + cos(-alpha)*cos(-gamma) )
  94. + cx*( -sin(-gamma)*cos(-beta) );
  95.  
  96. rotx=cz*(sin(-alpha)*sin(-gamma) - sin(-beta)*cos(-alpha)*cos(-gamma) )
  97. + cy*(sin(-gamma)*cos(-alpha) + sin(-alpha)*sin(-beta)*cos(-gamma) )
  98. + cx*(cos(-beta)*cos(-gamma));
  99. );
  100.  
  101. # Convert the rotated point BACK into spherical coordinates
  102. # Atan2 is a special programming thing to do arctan, except better
  103. # It's got its own Wikipedia page
  104.  
  105. ntheta=atan2(roty, rotx);
  106. nphi=acos( rotz / sqrt(rotx*rotx + roty*roty + rotz*rotz) );
  107.  
  108. # Convert spherical coordinates almost straight back into
  109. # equirectangular pixel coordinates.
  110.  
  111. x2=( (ntheta/pi/2)*w+w/2 );
  112. y2=(nphi/pi)*h;
  113.  
  114. # I: from image #0 take the vector(pixel value) of x2,y2,z(=0)
  115. # "2": Cubic interpolation of pixels
  116. # "0": Boundary condition: if you hit "nothing", what do you do? N/A here.
  117.  
  118. I(#0,x2,y2,0,2,1)"
  119.  
  120. # Keep the last image generated ("k." = keep[-1] image index)
  121.  
  122. -keep[-1]
  123.  
  124.  
  125. ########################
  126. ## SINUSOIDAL
  127. #@gmic sinusoidal_map : yaw angle, pitch angle, roll angle, bgred, bggreen, bgblue, slices
  128. #@gmic : Take a map or equirectangular panorama
  129. #@gmic : and rotate it around z, y and x angles,
  130. #@gmic : then project it as a sliced sinusoidal map
  131.  
  132. #@gui Sinusoidal Map : sinusoidal_map
  133. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a sinusoidal projection.")
  134. #@gui : Roll = float(0,-360,360)
  135. #@gui : Pitch ("north and south") = float(0,-360,360)
  136. #@gui : Yaw ("east and west") = float(0,-360,360)
  137. #@gui : Background color = color(128,128,128)
  138. #@gui : Background opacity (%) = float(100,0,100)
  139. #@gui : Slices = int(1,1,72)
  140. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  141.  
  142.  
  143. sinusoidal_map :
  144.  
  145. # $1-3 Default to 0 degrees
  146. # $4-6 color(r,g,b) gives you three separate parameters instead of, say, a vector
  147. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=0},${7=100} -check ${8=1}>=1
  148.  
  149. slices=$8
  150. # Give input alpha-channel to prevent errors with transparency
  151. to_a[0]
  152.  
  153. _fF_zoomaccurate=1
  154. _fF_zoomfactor=1
  155.  
  156. # First rotate the map with the appropriate custom command
  157. if $1!=0" || "$2!=0" || "$3!=0
  158. rotate_equirectangular_map $1,$2,$3
  159. fi
  160.  
  161. # Slice the rotated map image into vertical sections
  162. split x,$slices
  163.  
  164. echo "slices: "
  165. echo $slices
  166.  
  167. # Loop over each slice. $! = number of images in the image list
  168. repeat $!
  169.  
  170. # No keyword: create an image. The ""-block is the mathematical function
  171. # whose result/return value at the end determines the value of the pixel
  172. # x, y we are iterating over.
  173. -input 100%,100%,1,100%,"
  174.  
  175. # Background color
  176. bgr="{$4}";
  177. bgg="{$5}";
  178. bgb="{$6}";
  179. bgopacity="{$7}/100*255";
  180.  
  181. # Half-accidental magical formula. Half-sin is used to center the coordinates somehow??
  182. halfsin=sin( (y/h)*pi )*w / 2;
  183. x2=(x+halfsin-w/2)/max(0.0000000000000000000000000000000000000001, sin( y/h*pi));
  184. # y = y in this case, just use original y.
  185.  
  186. # Check whether we're coloring inside or outside the sinusoidal shape
  187. if ( (x >= (w/2-halfsin) && x <= (w/2+halfsin) ),
  188. I(#"$>",x2,y,0,2,0),
  189. color=[bgr, bgg, bgb, bgopacity]);"
  190. # repeat-loop end
  191. done
  192.  
  193. echo "AFTER REPEAT\n"
  194.  
  195. # Delete the old map, which was sliced up
  196. remove[0-{$slices-1}]
  197.  
  198. # join all the rest (new) slices together in order on the x axis
  199. append x
  200.  
  201. # Keep the last (. = [-1] image index) image generated
  202. -keep[-1]
  203.  
  204.  
  205.  
  206. #######################
  207. ## TRIANGULAR
  208. #@gmic triangular_projection : yaw angle, pitch angle, roll angle, bgred, bggreen, bgblue, slices
  209. #@gmic : Take a map or equirectangular panorama
  210. #@gmic : and rotate it around z, y and x angles,
  211. #@gmic : then project it as a sliced triangular map.
  212. #@gmic : This is not the Collignon projection.
  213. #@gmic : It's a naive triangular interpolation.
  214.  
  215. #@gui Triangular Projection : triangular_projection
  216. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a naive rectilinear projection. This is not a Collignon projection, but a naive linear interpolation.")
  217. #@gui : note = note("When the width at the poles is 50%, this becomes the Eckert I projection. These projections are not equal-area nor do they have any other useful cartographical property.")
  218. #@gui : Roll = float(0,-360,360)
  219. #@gui : Pitch ("north and south") = float(0,-360,360)
  220. #@gui : Yaw ("east and west") = float(0,-360,360)
  221. #@gui : Background color = color(128,128,128)
  222. #@gui : Background opacity (%) = float(100,0,100)
  223. #@gui : Slices = int(1,1,72)
  224. #@gui : Width at the poles % = float(50,0,100)
  225. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  226.  
  227.  
  228. triangular_projection :
  229.  
  230.  
  231. # $1-3 Default to 0 degrees
  232. # $4-7 color(r,g,b) gives you three separate parameters instead of, say, a vector
  233. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=0},${7=100} -check ${8=1}>=1" && "${9=5}>=0
  234.  
  235. slices=$8
  236. magic=$9
  237.  
  238. # Give input alpha-channel to prevent errors with transparency
  239. to_a[0]
  240.  
  241.  
  242. # Background color
  243. bgr=$4;
  244. bgg=$5;
  245. bgb=$6;
  246. bgopacity={$7/100*255};
  247.  
  248. # Give input alpha-channel to prevent errors with transparency
  249. to_a[0]
  250.  
  251. _fF_zoomaccurate=1
  252. _fF_zoomfactor=1
  253.  
  254. # First rotate the map with the appropriate custom command
  255. if $1!=0" || "$2!=0" || "$3!=0
  256. rotate_equirectangular_map $1,$2,$3
  257. fi
  258.  
  259. # Slice the rotated map image into vertical sections
  260. split x,$slices
  261.  
  262. # Input: create an image. The ""-block is the mathematical function
  263. # whose result/return value at the end determines the value of the pixel
  264. # x, y we are iterating over.
  265. -input 100%,100%,1,4,"
  266.  
  267.  
  268. startwidth="$magic";
  269.  
  270. if( (y<=h/2),
  271. my=y,
  272. my=h-y);
  273.  
  274. sc=startwidth/100;
  275. squish=w / lerp((w*sc), w, my/(h/2));
  276. shift=(w-squish*w)/2;
  277. X = squish * x +shift;
  278. Y=y;
  279.  
  280. #debug normalisation values: cap at w so everything normalizes
  281. #result=[min(w,X),min(Y,w),0,w];"
  282.  
  283. result=[X,Y,0,0];"
  284.  
  285. # Loop over each slice. $! = number of images in the image list
  286. # Because the warp map is identical for each "slice" in this version
  287. # of the script / process, it doesn't need to generated for each slice,
  288. # so it's just reused. The sinusoidal script does an in-place operation
  289. # for each slice.
  290.  
  291. repeat $slices
  292. # the $> is a loop index incrementing forwards (think i++).
  293. # if there was $< it would loop backwards (like i--)
  294. # but here we only need to use the one warp image at the very
  295. # end "[-1]" of the image stack
  296. -warp[$>] [-1],0,2,0
  297. done
  298.  
  299. #DISPLACEMENT WARP
  300. # Delete the warp map, which was sliced up
  301. remove[-1]
  302.  
  303. # join all the rest (new) slices together in order on the x axis
  304.  
  305. #for debugging by normalizing and appending the warp field image
  306. #-normalize[-1] 0,255
  307.  
  308. append x
  309.  
  310. # generate blank image, then fill_color with bg-color
  311. -input 100%,100%,1,4,0
  312. -fill_color[-1] $bgr,$bgg,$bgb,$bgopacity
  313.  
  314. reverse
  315. -blend alpha
  316.  
  317. # Keep the last (. = [-1] image index) image generated
  318. -keep[-1]
  319.  
  320.  
  321. #######################
  322. ## CYLINDRICAL EQUAL AREA
  323. #@gmic cylindrical_equal_area : yaw angle, pitch angle, roll angle, standard_parallel angle
  324. #@gmic : Take a map or equirectangular panorama
  325. #@gmic : and rotate it around z, y and x angles,
  326. #@gmic : then project it as a cylindrical equal area map.
  327.  
  328. #@gui Cylindrical Equal-Area Projection : cylindrical_equal_area
  329. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a cylindrical equal-area projection.")
  330. #@gui : Roll = float(0,-360,360)
  331. #@gui : Pitch ("north and south") = float(0,-360,360)
  332. #@gui : Yaw ("east and west") = float(0,-360,360)
  333. #@gui : Standard parallels = float(45,0,80)
  334. #@gui : Named specializations = choice("[Standard parallels ↑]","Lambert 0° π:1","Behrmann 30°","Smyth/Craster ≈37°04\′17\″ 2:1","Trystan Edwards 37°24\′","Hobo-Dyer 37°30\′","Gall-Peters 45°","Balthasar 50°","Tobler ≈55°39\′14\″ 1:1","[Custom formula ↓]")
  335. #@gui : Custom formula = text{"acos(sqrt(1/pi))"}
  336. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  337.  
  338.  
  339. cylindrical_equal_area :
  340.  
  341. # $1-3 Default to 0 degrees
  342. -skip ${1=0},${2=0},${3=0} -check ${4=45}<=80
  343.  
  344. # Give input alpha-channel to prevent errors with transparency
  345. to_a[0]
  346.  
  347. # Standard parallells custom
  348. phi_input=$4
  349.  
  350. # G'mic-style "array"
  351. if $5==1 phi0=0
  352. elif $5==2 phi0={30*pi/180}
  353. elif $5==3 phi0={acos(sqrt(2/pi))}
  354. elif $5==4 phi0={37.4*pi/180}
  355. elif $5==5 phi0={37.5*pi/180}
  356. elif $5==6 phi0={45*pi/180}
  357. elif $5==7 phi0={50*pi/180}
  358. elif $5==8 phi0={acos(sqrt(1/pi))}
  359. elif $5==9 phi0={$6}
  360. # use Standards parallels input
  361. else phi0={$phi_input*pi/180}
  362. fi
  363. # phi0 is the input parameter Standard parallels in radians
  364.  
  365.  
  366.  
  367.  
  368. _fF_zoomaccurate=1
  369. _fF_zoomfactor=1
  370.  
  371. # First rotate the map with the appropriate custom command
  372. if $1!=0" || "$2!=0" || "$3!=0
  373. rotate_equirectangular_map $1,$2,$3
  374. fi
  375.  
  376.  
  377. # width of original map
  378. nw={0,w}
  379.  
  380. # height of original map
  381. nh={0,h}
  382.  
  383. # Determine the dimensions of the output projection by using
  384. # the ratio of w:h as give by Wikipedia
  385. # https://en.wikipedia.org/wiki/Cylindrical_equal-area_projection
  386. # w:h = pi*(cos phi_0)²
  387. W={$nw*pi*cos($phi0)^2}
  388. H={$nw}
  389.  
  390. # Initialize the image using the above precalculated width and height
  391. # The "" block contains the mathematical 'formula' that determines
  392. # what each x,y pixel is.
  393.  
  394. -input $W,$H,1,4,"
  395.  
  396. phi0="$phi0";
  397. pi2=2*pi;
  398.  
  399. # Use centered coordinates. 0,0 is in the middle of the
  400. # projection
  401. cx=x-w/2;
  402. cy=y-h/2;
  403.  
  404. # Turn them into radians??
  405. cx = (cx/w *pi2);
  406. cy = (cy/w *pi2);
  407.  
  408. # Somehow I cancelled out the formula for nX, it's
  409. # now part of the output image width formula above????
  410. nX = cx;
  411. # I don't know! I don't know! But this formula works!!
  412. nY = asin(cy * cos(phi0)^2);
  413.  
  414. # Here's the magic trick: the X and Y coordinates we are
  415. # calculating are not about where something goes, it's about
  416. # where something COMES FROM. Because of the way G'mic makes
  417. # you think about 'filters', you have to invert the transformation
  418. # direction: you look at the OUTPUT images' x and y coordinates,
  419. # then figure out where *from* that values has to come from.
  420. # Basically, we're taking the output and calculating it in reverse,
  421. # so that we get equirectangular (=spherical) coordinates back.
  422. # That's why instead of calculating with SINE, which is the way you
  423. # do it forward, we calculate with ARCSINE, which is backwards.
  424. # This is the number one reason why doing this stuff in G'mic is torture.
  425. # Check out the formulas in
  426. # https://mathworld.wolfram.com/CylindricalEqual-AreaProjection.html
  427.  
  428. # EQUIRECTANGULAR COORDINATES
  429. # Get the width and height of the original equirectangular input
  430. ow="{0,w}";
  431. oh="{0,h}";
  432.  
  433. # Pretty simply translate radians coordinates straight into
  434. # equirectangular pixel coordinates.
  435. X=ow/2 + (ow * nX/pi2) ;
  436. Y=oh/2 + (oh * nY/pi) ;
  437.  
  438. if( (nX>=-pi && nX<=pi && nY>=-pi/2 && nY<=pi/2),
  439. # TAKE the value of X,Y FROM the original equirectangular source map
  440. # and input the value into x,y in the output.
  441. result=I(#0,X%ow,Y,0,2,0),
  442. # else output transparent pixel
  443. result=[0,0,0,0]
  444. );
  445.  
  446. result"
  447.  
  448. -autocrop[1]
  449.  
  450. # Keep the last (. = [-1] image index) image generated
  451. -keep[-1]
  452.  
  453.  
  454.  
  455. #######################
  456. ## LAMBERT AZIMUTHAL
  457. #@gmic lambert_azimuthal_projection : yaw angle, pitch angle, roll angle
  458. #@gmic : Take a map or equirectangular panorama
  459. #@gmic : and rotate it around z, y and x angles,
  460. #@gmic : then project it as a Lambert azimuthal equal area map.
  461.  
  462. #@gui Lambert Azimuthal Equal-Area : lambert_azimuthal_projection
  463. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a Lambert Azimuthal Equal-Area map.")
  464. #@gui : Roll = float(0,-360,360)
  465. #@gui : Pitch = float(0,-360,360)
  466. #@gui : Yaw = float(0,-360,360)
  467. #@gui : Central longitude = float(0,-180,180)
  468. #@gui : Standard parallel = float(0,-90,90)
  469. #@gui : Cut-off at = float(180,0,180)
  470. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  471.  
  472.  
  473. lambert_azimuthal_projection :
  474.  
  475. # $1-3 Default to 0 degrees
  476. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=180}
  477.  
  478. phi0=0*pi/180
  479. cutoff={$6*pi/180}
  480.  
  481. # Give input alpha-channel to prevent errors with transparency
  482. to_a[0]
  483.  
  484.  
  485. _fF_zoomaccurate=1
  486. _fF_zoomfactor=1
  487.  
  488. # First rotate the map with the appropriate custom command
  489. if $1!=0" || "$2!=0" || "$3!=0
  490. rotate_equirectangular_map $1,$2,$3
  491. fi
  492.  
  493. # There is a bug in the naive formulas I put into this
  494. # script. Changing phi0 causes the map to get a literal
  495. # rip hole into it. If that was the only way for me to
  496. # change the angle of the globe view, this script
  497. # would be basically pretty useless and would only make
  498. # maps in one orientation. Happily, you can subsitute
  499. # those calculations with the usual equirectangular
  500. # rotation command. In this case, the two below are
  501. # needed to make "parallel" and "longitude" make sense.
  502. # The equirectangular_projection command applies the
  503. # rotations in a specific order. If you try to use it
  504. # *like* longitudes and parallels, you have to find the
  505. # equivalent values somehow. Below, to easily get the
  506. # "right" orientation, first East-West movement is
  507. # applied, then North-South.
  508. if $4!=0
  509. rotate_equirectangular_map 0,0,$4
  510. fi
  511.  
  512. if $5!=0
  513. rotate_equirectangular_map 0,{-$5},0
  514. fi
  515.  
  516.  
  517. # width of original map
  518. nw={0,w}
  519.  
  520. # height of original map
  521. nh={0,h}
  522.  
  523. # Initialize the image using the above precalculated width and height
  524. # The "" block contains the mathematical 'formula' that determines
  525. # what each x,y pixel is.
  526.  
  527. -input {$nw+1},{$nw+1},1,4,"
  528.  
  529. #inside the ""-block it's easier to pre-translate dollar-variables
  530. phi0="$phi0";
  531. cutoff="$cutoff";
  532. pi2=2*pi;
  533.  
  534.  
  535. # Use centered coordinates. 0,0 is in the middle of the
  536. # projection
  537. cx=x-w/2;
  538. cy=y-h/2;
  539.  
  540.  
  541. # Turn them into radians??
  542. cx = (cx/w *pi2);
  543. cy = (cy/w *pi2);
  544.  
  545. # For convenience sake. Rho is obviously, in retrospect
  546. # distance (from center)... At rho==sqrt/2 is where the
  547. # switch in the atan-function should happen.
  548. # These are from
  549. # https://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html
  550. rho=sqrt(cx^2 + cy^2);
  551. C=2*asin(1/2*rho);
  552.  
  553. # atan2 can be used here! The formula looks complicated, but just
  554. # take the parts of the fraction and replace the / with , for the function
  555. nX = atan2( (cx*sin(C)), ( (rho*cos(phi0)*cos(C))-(cy*sin(phi0)*sin(C))) );
  556.  
  557. nY = asin( cos(C)*sin(phi0)+(cy*sin(C)*cos(phi0)) / rho);
  558. # the sin(phi0)'s cancel out with multiplying by zero, should remove...
  559.  
  560. # For the logic behind inverting the formula, check out the
  561. # Cylindrical Equal-Area filter for a rant
  562.  
  563. # EQUIRECTANGULAR COORDINATES
  564. # Get the width and height of the original equirectangular input
  565. ow="{0,w}";
  566. oh="{0,h}";
  567.  
  568. # Pretty simply translate radians coordinates straight into
  569. # equirectangular pixel coordinates.
  570. X=ow/2 + (ow * nX/pi2);
  571. Y=oh/2 + (oh * nY/pi);
  572.  
  573. # Calculate the map's rendering radius based on user's cutoff angle
  574. # This is a simplified Great Circle distance formula
  575. # where half the formula disappears because sin(0)=0
  576. deltaS=acos(cos(nY)*cos(nX));
  577.  
  578. # check whether the pixel's geographical nX, nY coordinates
  579. # are within the normal range AND inside the deltaS cutoff
  580. if( (nX>=-pi && nX<=pi && nY>=-pi/2 && nY<=pi/2 && abs(deltaS)<=cutoff),
  581. # TAKE the value of X,Y FROM the original equirectangular source map
  582. # and input the value into x,y in the output.
  583. result=I(#0,X,Y,0,2,0),
  584. # else output transparent pixel
  585. result=[0,0,0,0]
  586. );
  587. # return the last value
  588. result;"
  589.  
  590. # New trick: if you can't calculate the width and height of the image
  591. # in advance, doing it on a "blank" square image can work as well.
  592. # In the ""-block of the -input above the last pixel assignment
  593. # checks if the geographical coordinates we're looking for make
  594. # any sense, that is, no weird values above or below plus-minus 180 or
  595. # 90 degrees. If the pixel x,y coordinates would result in something
  596. # nonsensical, output a transparent pixel. This is necessary, because it's
  597. # not automatic: you get weird blocks of white *and* transparent
  598. # lines. Then you can use autocrop from g'mic to do some magic and remove
  599. # the transparent areas.
  600. -autocrop[-1]
  601.  
  602.  
  603. # Keep the last (. = [-1] image index) image generated
  604. -keep[-1]
  605.  
  606.  
  607. #######################
  608. ## MOLLWEIDE
  609. #@gmic mollweide_projection : yaw angle, pitch angle, roll angle, central_longitude angle, standard_parallel angle
  610. #@gmic : Take a map or equirectangular panorama
  611. #@gmic : and rotate it around z, y and x angles,
  612. #@gmic : then project it as a Mollweide projection map.
  613.  
  614. #@gui Mollweide : mollweide_projection
  615. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a Mollweide equal-area pseudocylindrical map.")
  616. #@gui : Roll = float(0,-360,360)
  617. #@gui : Pitch = float(0,-360,360)
  618. #@gui : Yaw = float(0,-360,360)
  619. #@gui : Central longitude = float(0,-180,180)
  620. #@gui : Standard parallel = float(0,-90,90)
  621. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  622.  
  623.  
  624. mollweide_projection :
  625.  
  626. # $1-3 Default to 0 degrees
  627. -skip ${1=0},${2=0},${3=0},${4=0},${5=0}
  628.  
  629. # Give input alpha-channel to prevent errors with transparency
  630. to_a[0]
  631.  
  632. _fF_zoomaccurate=1
  633. _fF_zoomfactor=1
  634.  
  635. # First rotate the map with the appropriate custom command
  636. if $1!=0" || "$2!=0" || "$3!=0
  637. rotate_equirectangular_map $1,$2,$3
  638. fi
  639.  
  640. # To apply more standard "longitude" and "standard parallel
  641. # translations, use the same rotate command, except in an order
  642. # where longitude (yaw) is applied first and latitude (pitch)
  643. # is applied second, in two different operations.
  644. if $4!=0
  645. rotate_equirectangular_map 0,0,$4
  646. fi
  647.  
  648. if $5!=0
  649. rotate_equirectangular_map 0,{-$5},0
  650. fi
  651.  
  652.  
  653. # width and height of original map
  654. nw={0,w}
  655. nh={0,h}
  656.  
  657. # Initialize the image using the above precalculated width and height
  658. # The "" block contains the mathematical 'formula' that determines
  659. # what each x,y pixel is.
  660.  
  661. -input {$nw+1},{$nw+1},1,4,"
  662.  
  663. #inside the ""-block it's easier to pre-translate dollar-variables
  664. pi2=2*pi;
  665.  
  666.  
  667. # Use centered coordinates. 0,0 is in the middle of the
  668. # projection
  669. cx=x-w/2;
  670. cy=y-h/2;
  671.  
  672.  
  673. # Turn them into radians
  674. cx = (cx/w *pi2);
  675. cy = (cy/w *pi2);
  676.  
  677.  
  678. # constants for the formulas. R=1 is the Radius of the globe,
  679. # I've put 1 here just because i couldn't find a convenient
  680. # fraction of pi, and it doesn't change the shape, just scale.
  681. # Theta is separated as a constant following the example of
  682. # the mathworld equations I found:
  683. # https://mathworld.wolfram.com/MollweideProjection.html
  684. R=1;
  685. theta=asin(cy/(sqrt(2)*R));
  686.  
  687. nX = (pi*cx) / (2*R*sqrt(2)*cos(theta));
  688. nY = asin( ((2*theta)+sin(2*theta)) / pi);
  689.  
  690.  
  691. # For the logic why I'm inverting the formula, check out the
  692. # Cylindrical Equal-Area filter for a rant
  693.  
  694. # EQUIRECTANGULAR COORDINATES
  695. # Get the width and height of the original equirectangular input
  696. ow="{0,w}";
  697. oh="{0,h}";
  698.  
  699. # Pretty simply translate radians coordinates straight into
  700. # equirectangular pixel coordinates.
  701. X=ow/2 + (ow * nX/pi2);
  702. Y=oh/2 + (oh * nY/pi);
  703.  
  704. # check whether the pixel's geographical nX, nY coordinates
  705. # are within the normal range
  706. if( (nX>=-pi && nX<=pi && nY>=-pi/2 && nY<=pi/2),
  707. # then TAKE the value of X,Y FROM the original equirectangular
  708. # source map and input the value into x,y in the output.
  709. result=I(#0,X,Y,0,2,0),
  710. # else output a transparent pixel
  711. result=[0,0,0,0]
  712. );
  713. # return the last value
  714. result;"
  715.  
  716. # New trick: if you can't calculate the width and height of the image
  717. # in advance, doing it on a "blank" square image can work as well.
  718. # In the ""-block of the -input above the last pixel assignment
  719. # checks if the geographical coordinates we're looking for make
  720. # any sense, that is, no weird values above or below plus-minus 180 or
  721. # 90 degrees. If the pixel x,y coordinates would result in something
  722. # nonsensical, output a transparent pixel. This is necessary, because it's
  723. # not automatic: you get weird blocks of white *and* transparent
  724. # lines. Then you can use autocrop from g'mic to do some magic and remove
  725. # the transparent areas.
  726. -autocrop[-1]
  727.  
  728.  
  729. # Keep the last (. = [-1] image index) image generated
  730. -keep[-1]
  731.  
  732.  
  733. #######################
  734. ## ORTHOGRAPHIC
  735. #@gmic orthographic_projection : yaw angle, pitch angle, roll angle, central_longitude angle, standard_parallel angle
  736. #@gmic : Take a map or equirectangular panorama
  737. #@gmic : and rotate it around z, y and x angles,
  738. #@gmic : then project it as an orthographic projection map.
  739.  
  740. #@gui Orthographic : orthographic_projection
  741. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as an orthographic hemisphere map. This map projection is basically the view you'd get if you looked at earth from an infinite distance.")
  742. #@gui : Roll = float(0,-360,360)
  743. #@gui : Pitch = float(0,-360,360)
  744. #@gui : Yaw = float(0,-360,360)
  745. #@gui : Central longitude = float(0,-180,180)
  746. #@gui : Standard parallel = float(0,-90,90)
  747. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  748.  
  749.  
  750. orthographic_projection :
  751.  
  752. # $1-3 Default to 0 degrees
  753. -skip ${1=0},${2=0},${3=0},${4=0},${5=0}
  754.  
  755. # Give input alpha-channel to prevent errors with transparency
  756. to_a[0]
  757.  
  758. _fF_zoomaccurate=1
  759. _fF_zoomfactor=1
  760.  
  761. # First rotate the map with the appropriate custom command
  762. if $1!=0" || "$2!=0" || "$3!=0
  763. rotate_equirectangular_map $1,$2,$3
  764. fi
  765.  
  766. # To apply more standard "longitude" and "standard parallel
  767. # translations, use the same rotate command, except in an order
  768. # where longitude (yaw) is applied first and latitude (pitch)
  769. # is applied second, in two different operations.
  770. if $4!=0
  771. rotate_equirectangular_map 0,0,$4
  772. fi
  773.  
  774. if $5!=0
  775. rotate_equirectangular_map 0,{-$5},0
  776. fi
  777.  
  778.  
  779. # width and height of original map
  780. nw={0,w}
  781. nh={0,h}
  782.  
  783. # Initialize the image using the above precalculated width and height
  784. # The "" block contains the mathematical 'formula' that determines
  785. # what each x,y pixel is.
  786.  
  787. -input {$nw+1},{$nw+1},1,4,"
  788.  
  789. #inside the ""-block it's easier to pre-translate dollar-variables
  790. pi2=2*pi;
  791.  
  792.  
  793. # Use centered coordinates. 0,0 is in the middle of the
  794. # projection
  795. cx=x-w/2;
  796. cy=y-h/2;
  797.  
  798.  
  799. # Turn them into radians
  800. cx = (cx/w *pi2);
  801. cy = (cy/w *pi2);
  802.  
  803.  
  804. # Formulas and constants:
  805. # https://en.wikipedia.org/wiki/Orthographic_projection_in_cartography
  806. # R, radius of globe seems to work best with R=pi/2
  807. R=pi/2;
  808. rho=sqrt(cx^2 + cy^2);
  809. C=asin(rho/R);
  810. phi0=0;
  811. l0=0;
  812.  
  813. # the formulas are simplified a bit because this script uses
  814. # phi0=0 and l0=0, leading to anything with sin(phi0) to
  815. # multiply with 0. The standard parallel and central longitude
  816. # are calculated earlier with actual global rotations.
  817. # Happily, here we can use atan2.
  818. nY = asin((cy*sin(C)*cos(phi0))/rho);
  819. nX = atan2((cx*sin(C)), (rho*cos(C)*cos(phi0)));
  820.  
  821. # For the logic why I'm inverting the formula, check out the
  822. # Cylindrical Equal-Area filter for a rant
  823.  
  824. # EQUIRECTANGULAR COORDINATES
  825. # Get the width and height of the original equirectangular input
  826. ow="{0,w}";
  827. oh="{0,h}";
  828.  
  829. # Pretty simply translate radians coordinates straight into
  830. # equirectangular pixel coordinates.
  831. X=ow/2 + (ow * nX/pi2);
  832. Y=oh/2 + (oh * nY/pi);
  833.  
  834. # check whether the pixel's geographical nX, nY coordinates
  835. # are within the normal range
  836. # Why does this break if the condition is set to true manually???
  837. if( (nX>-pi && nX<pi && nY>-pi/2 && nY<pi/2),
  838. # then TAKE the value of X,Y FROM the original equirectangular
  839. # source map and input the value into x,y in the output.
  840. result=I(#0,X,Y,0,2,0),
  841. # else output a transparent pixel
  842. result=[0,0,0,0]
  843. );
  844. # return the last value
  845. result;"
  846.  
  847. # New trick: if you can't calculate the width and height of the image
  848. # in advance, doing it on a "blank" square image can work as well.
  849. # In the ""-block of the -input above the last pixel assignment
  850. # checks if the geographical coordinates we're looking for make
  851. # any sense, that is, no weird values above or below plus-minus 180 or
  852. # 90 degrees. If the pixel x,y coordinates would result in something
  853. # nonsensical, output a transparent pixel. This is necessary, because it's
  854. # not automatic: you get weird blocks of white *and* transparent
  855. # lines. Then you can use autocrop from g'mic to do some magic and remove
  856. # the transparent areas.
  857. -autocrop[-1]
  858.  
  859.  
  860. # Keep the last (. = [-1] image index) image generated
  861. -keep[-1]
  862.  
  863. #######################
  864. ## AZIMUTHAL EQUIDISTANT
  865. #@gmic azimuthal_equidistant_projection : yaw angle, pitch angle, roll angle, central_longitude angle, standard_parallel angle
  866. #@gmic : Take a map or equirectangular panorama
  867. #@gmic : and rotate it around z, y and x angles,
  868. #@gmic : then project it as an azimuthal equidistant map.
  869.  
  870. #@gui Azimuthal Equidistant : azimuthal_equidistant_projection
  871. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as an azimuthal equidistant map. All points on the map are at proportionally correct distances from the center point, and that all points on the map are at the correct azimuth (direction) from the center point.")
  872. #@gui : Roll = float(0,-360,360)
  873. #@gui : Pitch = float(0,-360,360)
  874. #@gui : Yaw = float(0,-360,360)
  875. #@gui : Central longitude = float(0,-180,180)
  876. #@gui : Standard parallel = float(0,-90,90)
  877. #@gui : Cut-off at = float(180,0,180)
  878. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  879.  
  880.  
  881. azimuthal_equidistant_projection :
  882.  
  883. # $1-3 Default to 0 degrees
  884. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=180}
  885.  
  886. # Give input alpha-channel to prevent errors with transparency
  887. to_a[0]
  888.  
  889. cutoff={$6*pi/180}
  890.  
  891. _fF_zoomaccurate=1
  892. _fF_zoomfactor=1
  893.  
  894. # First rotate the map with the appropriate custom command
  895. if $1!=0" || "$2!=0" || "$3!=0
  896. rotate_equirectangular_map $1,$2,$3
  897. fi
  898.  
  899. # To apply more standard "longitude" and "standard parallel
  900. # translations, use the same rotate command, except in an order
  901. # where longitude (yaw) is applied first and latitude (pitch)
  902. # is applied second, in two different operations.
  903. if $4!=0
  904. rotate_equirectangular_map 0,0,$4
  905. fi
  906.  
  907. if $5!=0
  908. rotate_equirectangular_map 0,{-$5},0
  909. fi
  910.  
  911.  
  912. # width and height of original map
  913. nw={0,w}
  914. nh={0,h}
  915.  
  916. # Initialize the image using the above precalculated width and height
  917. # The "" block contains the mathematical 'formula' that determines
  918. # what each x,y pixel is.
  919.  
  920. -input {$nw+1},{$nw+1},1,4,"
  921.  
  922. #inside the ""-block it's easier to pre-translate dollar-variables
  923. pi2=2*pi;
  924. cutoff="$cutoff";
  925.  
  926.  
  927. # Use centered coordinates. 0,0 is in the middle of the
  928. # projection
  929. cx=x-w/2+1;
  930. cy=y-h/2+1;
  931.  
  932.  
  933. # Turn them into radians
  934. cx = (cx/w *pi2);
  935. cy = (cy/w *pi2);
  936.  
  937.  
  938. # Formulas and constants:
  939. # https://en.wikipedia.org/wiki/Orthographic_projection_in_cartography
  940. # R, radius of globe seems to work best with R=pi/2
  941. C=sqrt(cx^2 + cy^2);
  942. phi0=0;
  943. l0=0;
  944.  
  945. # the formulas are simplified a bit because this script uses
  946. # phi0=0 and l0=0, leading to anything with sin(phi0) to
  947. # multiply with 0. The standard parallel and central longitude
  948. # are calculated earlier with actual global rotations.
  949. # Happily, here we can use atan2.
  950. nY = asin( (cy*sin(C)) / C );
  951. nX = atan2((cx*sin(C)), (C*cos(C)));
  952.  
  953. # For the logic why I'm inverting the formula, check out the
  954. # Cylindrical Equal-Area filter for a rant
  955.  
  956. # EQUIRECTANGULAR COORDINATES
  957. # Get the width and height of the original equirectangular input
  958. ow="{0,w}";
  959. oh="{0,h}";
  960.  
  961. # Pretty simply translate radians coordinates straight into
  962. # equirectangular pixel coordinates.
  963. X=ow/2 + (ow * nX/pi2);
  964. Y=oh/2 + (oh * nY/pi);
  965.  
  966. # Calculate the map's rendering radius based on user's cutoff angle
  967. # This is a simplified Great Circle distance formula
  968. # where half the formula disappears because sin(0)=0
  969. deltaS=acos(cos(nY)*cos(nX));
  970.  
  971. # check whether the pixel's geographical nX, nY coordinates
  972. # are within the normal range AND inside the deltaS cutoff
  973. if( (nX>=-pi && nX<=pi && nY>=-pi/2 && nY<=pi/2 && abs(deltaS)<=cutoff),
  974. # TAKE the value of X,Y FROM the original equirectangular source map
  975. # and input the value into x,y in the output.
  976. result=I(#0,X,Y,0,2,0),
  977. # else output transparent pixel
  978. result=[0,0,0,0]
  979. );
  980. result"
  981.  
  982.  
  983. # generate blank image, then fill_color with bg-color
  984. -input 100%,100%,1,4,0
  985. -fill_color[-1] 0,0,0,0
  986. -shape_circle {-1,w}
  987. -normalize[-1] 0,255
  988. -negate[-1]
  989. +blend_fade[-2,-3] [-1]
  990.  
  991.  
  992. # New trick: if you can't calculate the width and height of the image
  993. # in advance, doing it on a "blank" square image can work as well.
  994. # In the ""-block of the -input above the last pixel assignment
  995. # checks if the geographical coordinates we're looking for make
  996. # any sense, that is, no weird values above or below plus-minus 180 or
  997. # 90 degrees. If the pixel x,y coordinates would result in something
  998. # nonsensical, output a transparent pixel. This is necessary, because it's
  999. # not automatic: you get weird blocks of white *and* transparent
  1000. # lines. Then you can use autocrop from g'mic to do some magic and remove
  1001. # the transparent areas.
  1002. -autocrop[-1]
  1003.  
  1004.  
  1005. # Keep the last (. = [-1] image index) image generated
  1006. -keep[-1]
  1007.  
  1008.  
  1009.  
  1010.  
  1011.  
  1012.  
  1013.  
  1014. #######################
  1015. ## CONIC EQUIDISTANT
  1016. #@gmic conic_equidistant_projection : yaw angle, pitch angle, roll angle, central_longitude angle, reference_parallel angle, standard_parallel1 angle, standard_parallel2 angle
  1017. #@gmic : Take a map or equirectangular panorama
  1018. #@gmic : and rotate it around z, y and x angles,
  1019. #@gmic : then project it as a conic equidistant map.
  1020.  
  1021. #@gui Conic Equidistant : conic_equidistant_projection
  1022. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a conic equidistant map. Neither equal-area or conformal.")
  1023. #@gui : Roll = float(0,-360,360)
  1024. #@gui : Pitch = float(0,-360,360)
  1025. #@gui : Yaw = float(0,-360,360)
  1026. #@gui : Central longitude = float(0,-180,180)
  1027. #@gui : Reference parallel = float(0,-90,90)
  1028. #@gui : Standard parallel #1 = float(20,-90,90)
  1029. #@gui : Standard parallel #2 = float(60,-90,90)
  1030. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  1031.  
  1032.  
  1033. conic_equidistant_projection :
  1034.  
  1035. # $1-3 Default to 0 degrees
  1036. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=20},${7=60}
  1037.  
  1038. cutoff={$6*pi/180}
  1039. l0=$4
  1040. phi0=$5
  1041. phi1=$6
  1042. phi2=$7
  1043.  
  1044. # Give input alpha-channel to prevent errors with transparency
  1045. to_a[0]
  1046.  
  1047. _fF_zoomaccurate=1
  1048. _fF_zoomfactor=1
  1049.  
  1050. # First rotate the map with the appropriate custom command
  1051. if $1!=0" || "$2!=0" || "$3!=0
  1052. rotate_equirectangular_map $1,$2,$3
  1053. fi
  1054.  
  1055. # To apply more standard "longitude" and "standard parallel
  1056. # translations, use the same rotate command, except in an order
  1057. # where longitude (yaw) is applied first and latitude (pitch)
  1058. # is applied second, in two different operations.
  1059. # if $4!=0
  1060. # rotate_equirectangular_map 0,0,$4
  1061. # fi
  1062. #
  1063. # if $5!=0
  1064. # rotate_equirectangular_map 0,{-$5},0
  1065. # fi
  1066.  
  1067.  
  1068. # width and height of original map
  1069. nw={0,w}
  1070. nh={0,h}
  1071.  
  1072. # Initialize the image using the above precalculated width and height
  1073. # The "" block contains the mathematical 'formula' that determines
  1074. # what each x,y pixel is.
  1075.  
  1076. -input {$nw+1},{$nw+1},1,4,"
  1077.  
  1078. #inside the ""-block it's easier to pre-translate dollar-variables
  1079. pi2=2*pi;
  1080. phi0="$phi0"*pi/180;
  1081. l0="$l0"*pi/180;
  1082.  
  1083. phi1="$phi1"*pi/180;
  1084. phi2="$phi2"*pi/180;
  1085.  
  1086.  
  1087.  
  1088.  
  1089. # Use centered coordinates. 0,0 is in the middle of the
  1090. # projection
  1091. cx=x-w/2+1;
  1092. cy=y-h/2+1;
  1093.  
  1094.  
  1095. # Turn them into radians
  1096. cx = (cx/w *pi2*4/3);
  1097. cy = -(cy/w *pi2*4/3);
  1098.  
  1099.  
  1100. # Formulas and constants:
  1101. # https://mathworld.wolfram.com/ConicEquidistantProjection.html
  1102.  
  1103. if( (phi1==phi2),
  1104. n=sin(phi1),
  1105. n=(cos(phi1)-cos(phi2)) / (phi2-phi1)
  1106. );
  1107. G=cos(phi1)/n + phi1;
  1108. rho0=G-phi0;
  1109.  
  1110. rho=sign(n)*sqrt( cx^2 + (rho0-cy)^2);
  1111.  
  1112. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1113. # page 114, if using ATAN2, the signs must be reversed.
  1114. # I KNEW THERE WAS SOMETHING WRONG WITH THE FORMULA
  1115. # I had to dig deep.
  1116. if( (n>0),
  1117. (theta=atan2(cx, (rho0-cy))),
  1118. (theta=atan2(-cx, (-rho0+cy)))
  1119. );
  1120.  
  1121. nY = G-rho;
  1122. nX = l0 + theta/n;
  1123.  
  1124. # For the logic why I'm inverting the formula, check out the
  1125. # Cylindrical Equal-Area filter for a rant
  1126.  
  1127. # EQUIRECTANGULAR COORDINATES
  1128. # Get the width and height of the original equirectangular input
  1129. ow="{0,w}";
  1130. oh="{0,h}";
  1131.  
  1132. # Pretty simply translate radians coordinates straight into
  1133. # equirectangular pixel coordinates.
  1134. X=ow/2 + (ow * nX/pi2);
  1135. Y=oh/2 + (oh * -nY/pi);
  1136.  
  1137.  
  1138. # check whether the pixel's geographical nX, nY coordinates
  1139. # are within the normal range, here adjusted for l0 and phi0
  1140. # which are kind of superfluous here, they don't do much anything.
  1141.  
  1142. L=nX-l0;
  1143. P=nY-phi0;
  1144. if( (L>=-pi && L<=pi && P>=-pi/2 && P<=pi/2),
  1145. # TAKE the value of X,Y FROM the original equirectangular source map
  1146. # and input the value into x,y in the output.
  1147. result=I(#0,X%ow,Y,0,2,0),
  1148. # else output transparent pixel
  1149. result=[0,0,0,0]
  1150. );
  1151.  
  1152. result"
  1153.  
  1154.  
  1155. # New trick: if you can't calculate the width and height of the image
  1156. # in advance, doing it on a "blank" square image can work as well.
  1157. # In the ""-block of the -input above the last pixel assignment
  1158. # checks if the geographical coordinates we're looking for make
  1159. # any sense, that is, no weird values above or below plus-minus 180 or
  1160. # 90 degrees. If the pixel x,y coordinates would result in something
  1161. # nonsensical, output a transparent pixel. This is necessary, because it's
  1162. # not automatic: you get weird blocks of white *and* transparent
  1163. # lines. Then you can use autocrop from g'mic to do some magic and remove
  1164. # the transparent areas.
  1165. autocrop[-1]
  1166.  
  1167.  
  1168. # Keep the last (. = [-1] image index) image generated
  1169. -keep[-1]
  1170.  
  1171.  
  1172.  
  1173. #######################
  1174. ## LAMBERT CONFORMAL CONIC
  1175. #@gmic lambert_conformal_conic_projection : yaw angle, pitch angle, roll angle, central_longitude angle, reference_parallel angle, standard_parallel1 angle, standard_parallel2 angle, latcutoff angle, longcutoff angle
  1176. #@gmic : Take a map or equirectangular panorama
  1177. #@gmic : and rotate it around z, y and x angles,
  1178. #@gmic : then project it as a Lambert conic conformal map.
  1179.  
  1180. #@gui Lambert Conformal Conic : lambert_conformal_conic_projection
  1181. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a Lambert conic conformal map.")
  1182. #@gui : Roll = float(0,-360,360)
  1183. #@gui : Pitch = float(0,-360,360)
  1184. #@gui : Yaw = float(0,-360,360)
  1185. #@gui : Central longitude = float(0,-180,180)
  1186. #@gui : Reference parallel = float(0,-90,90)
  1187. #@gui : Standard parallel #1 = float(20,-90,90)
  1188. #@gui : Standard parallel #2 = float(60,-90,90)
  1189. #@gui : Crop latitude to = float(180,0,180)
  1190. #@gui : Crop longitude to = float(360,0,360)
  1191. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  1192.  
  1193.  
  1194. lambert_conformal_conic_projection :
  1195.  
  1196. # $1-3 Default to 0 degrees
  1197. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=20},${7=60},${8=180},${9=360}
  1198.  
  1199. latcutoff={$8*pi/180}
  1200. longcutoff={$9*pi/180}
  1201.  
  1202. l0=$4
  1203. phi0=$5
  1204. phi1=$6
  1205. phi2=$7
  1206.  
  1207. # Give input alpha-channel to prevent errors with transparency
  1208. to_a[0]
  1209.  
  1210. _fF_zoomaccurate=1
  1211. _fF_zoomfactor=1
  1212.  
  1213. # First rotate the map with the appropriate custom command
  1214. if $1!=0" || "$2!=0" || "$3!=0
  1215. rotate_equirectangular_map $1,$2,$3
  1216. fi
  1217.  
  1218. # width and height of original map
  1219. nw={0,w}
  1220. nh={0,h}
  1221.  
  1222. # Initialize the image using the above precalculated width and height
  1223. # The "" block contains the mathematical 'formula' that determines
  1224. # what each x,y pixel is.
  1225.  
  1226. -input {$nw+1},{$nw+1},1,4,"
  1227.  
  1228. #inside the ""-block it's easier to pre-translate dollar-variables
  1229. pi2=2*pi;
  1230. phi0="$phi0"*pi/180;
  1231. l0="$l0"*pi/180;
  1232.  
  1233. phi1="$phi1"*pi/180;
  1234. phi2="$phi2"*pi/180;
  1235.  
  1236. latcutoff="$latcutoff";
  1237. longcutoff="$longcutoff";
  1238.  
  1239.  
  1240.  
  1241.  
  1242. # Use centered coordinates. 0,0 is in the middle of the
  1243. # projection
  1244. cx=x-w/2+1;
  1245. cy=y-h/2+1;
  1246.  
  1247.  
  1248. # Turn them into radians
  1249. cx = (cx/w *pi2*pi/2);
  1250. cy = -(cy/w *pi2*pi/2);
  1251.  
  1252.  
  1253. # Formulas and constants:
  1254. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1255.  
  1256. if((phi1==phi2),
  1257. n=sin(phi1),
  1258. n= log(cos(phi1)/cos(phi2)) / log(tan((1/4)*pi+(1/2)*phi2)/tan((1/4)*pi+(1/2)*phi1))
  1259. );
  1260.  
  1261. F=(cos(phi1)*tan((1/4)*pi+(1/2)*phi1)^n) / n;
  1262.  
  1263. rho0=F/tan(1/4*pi+1/2*phi0)^n;
  1264.  
  1265. rho=sign(n)*sqrt( cx^2 + (rho0-cy)^2);
  1266.  
  1267. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1268. # page 114, if using ATAN2, the signs must be reversed.
  1269.  
  1270. if( (n>0),
  1271. (theta=atan2(cx, (rho0-cy))),
  1272. (theta=atan2(-cx, (-rho0+cy)))
  1273. );
  1274.  
  1275. # When the reference parallel is in the southern hemisphere,
  1276. # flip the projection and map
  1277. if( (phi0<0),
  1278. (n=-n)
  1279. );
  1280.  
  1281.  
  1282. if( (rho==0),
  1283. (nY = sign(n)*90*pi/180),
  1284. (nY = 2*atan((F/rho)^(1/n))-1/2*pi)
  1285. );
  1286. nX = l0 + theta/n;
  1287.  
  1288. # For the logic why I'm using inverse formulas, check out the
  1289. # Cylindrical Equal-Area filter
  1290.  
  1291. # EQUIRECTANGULAR COORDINATES
  1292. # Get the width and height of the original equirectangular input
  1293. ow="{0,w}";
  1294. oh="{0,h}";
  1295.  
  1296. # Pretty simply translate radians coordinates straight into
  1297. # equirectangular pixel coordinates.
  1298. X=ow/2 + (ow * nX/pi2);
  1299. Y=oh/2 + (oh * -nY/pi);
  1300.  
  1301.  
  1302. # check whether the pixel's geographical nX, nY coordinates
  1303. # are within the normal range, then further crop to the angles
  1304. # given by the crop values.
  1305.  
  1306. L=nX;
  1307. P=nY;
  1308. if( (L>=l0-longcutoff/2 && L<=l0+longcutoff/2 && P>=phi0-latcutoff/2 && P<=phi0+latcutoff/2),
  1309. # TAKE the value of X,Y FROM the original equirectangular source map
  1310. # and input the value into x,y in the output.
  1311. result=I(#0,X%ow,Y,0,2,0),
  1312. # else output transparent pixel
  1313. result=[0,0,0,0]
  1314. );
  1315.  
  1316. result"
  1317.  
  1318. # If the reference parallel is in the southern hemisphere,
  1319. # we've flipped the projection. Here we need to further rotate
  1320. # the image so it's not upside down.
  1321. if $phi0<0
  1322. rotate[-1] 180
  1323. fi
  1324.  
  1325.  
  1326. # New trick: if you can't calculate the width and height of the image
  1327. # in advance, doing it on a "blank" square image can work as well.
  1328. # In the ""-block of the -input above the last pixel assignment
  1329. # checks if the geographical coordinates we're looking for make
  1330. # any sense, that is, no weird values above or below plus-minus 180 or
  1331. # 90 degrees. If the pixel x,y coordinates would result in something
  1332. # nonsensical, output a transparent pixel. This is necessary, because it's
  1333. # not automatic: you get weird blocks of white *and* transparent
  1334. # lines. Then you can use autocrop from g'mic to do some magic and remove
  1335. # the transparent areas.
  1336. autocrop[-1]
  1337.  
  1338.  
  1339. # Keep the last (. = [-1] image index) image generated
  1340. -keep[-1]
  1341.  
  1342.  
  1343. #######################
  1344. ## MERCATOR
  1345. #@gmic mercator : yaw angle, pitch angle, roll angle, central_longitude angle, reference_parallel angle, standard_parallel1 angle, standard_parallel2 angle, latcutoff angle, longcutoff angle
  1346. #@gmic : Take a map or equirectangular panorama
  1347. #@gmic : and rotate it around z, y and x angles,
  1348. #@gmic : then project it as a Mercator projection.
  1349. #@gui Mercator Projection : mercator
  1350. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a Mercator projection. You can get a Transverse Mercator by rolling the globe by 90")
  1351. #@gui : Roll = float(0,-360,360)
  1352. #@gui : Pitch = float(0,-360,360)
  1353. #@gui : Yaw = float(0,-360,360)
  1354. #@gui : Central longitude = float(0,-180,180)
  1355. #@gui : Crop latitude to = float(85,0,85)
  1356. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  1357.  
  1358.  
  1359. mercator :
  1360.  
  1361. # $1-3 Default to 0 degrees
  1362. -skip ${1=0},${2=0},${3=0},${4=0},${5=85}
  1363.  
  1364. latcutoff={$5*pi/180}
  1365.  
  1366. l0=$4
  1367.  
  1368. # Give input alpha-channel to prevent errors with transparency
  1369. to_a[0]
  1370.  
  1371. _fF_zoomaccurate=1
  1372. _fF_zoomfactor=1
  1373.  
  1374. # First rotate the map with the appropriate custom command
  1375. if $1!=0" || "$2!=0" || "$3!=0
  1376. rotate_equirectangular_map $1,$2,$3
  1377. fi
  1378.  
  1379. # width and height of original map
  1380. nw={0,w}
  1381. nh={0,h}
  1382.  
  1383. # Initialize the image using the above precalculated width and height
  1384. # The "" block contains the mathematical 'formula' that determines
  1385. # what each x,y pixel is.
  1386.  
  1387. -input {$nw+1},{$nw+1},1,4,"
  1388.  
  1389. #inside the ""-block it's easier to pre-translate dollar-variables
  1390. pi2=2*pi;
  1391. l0="$l0"*pi/180;
  1392.  
  1393. latcutoff="$latcutoff";
  1394.  
  1395. # Use centered coordinates. 0,0 is in the middle of the
  1396. # projection
  1397. cx=x-w/2+1;
  1398. cy=y-h/2+1;
  1399.  
  1400.  
  1401. # Turn them into radians
  1402. cx = (cx/w *pi2);
  1403. cy = -( (cy/w *pi2) );
  1404.  
  1405.  
  1406. # Formulas and constants:
  1407. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1408.  
  1409. nY=atan(sinh(cy));
  1410. nX=l0+cx;
  1411.  
  1412. # For the logic why I'm using inverse formulas, check out the
  1413. # Cylindrical Equal-Area filter
  1414.  
  1415. # EQUIRECTANGULAR COORDINATES
  1416. # Get the width and height of the original equirectangular input
  1417. ow="{0,w}";
  1418. oh="{0,h}";
  1419.  
  1420. # Pretty simply translate radians coordinates straight into
  1421. # equirectangular pixel coordinates.
  1422. X=ow/2 + (ow * nX/pi2);
  1423. Y=oh/2 + (oh * -nY/pi);
  1424.  
  1425.  
  1426. # check whether the pixel's geographical nX, nY coordinates
  1427. # are within the normal range, then further crop to the angles
  1428. # given by the crop values.
  1429.  
  1430. L=nX-l0;
  1431. P=nY;
  1432. if( (L>=-pi && L<=pi && P>=-latcutoff && P<=+latcutoff),
  1433. # TAKE the value of X,Y FROM the original equirectangular source map
  1434. # and input the value into x,y in the output.
  1435. # X%ow, when combined with an l0 longitude shift, allows for
  1436. # rotating around the globe, otherwise it just cuts off at 180.
  1437. result=I(#0,X%ow,Y,0,2,0),
  1438. # else output transparent pixel
  1439. result=[0,0,0,0]
  1440. );
  1441.  
  1442. result"
  1443.  
  1444. autocrop[-1]
  1445.  
  1446.  
  1447. # Keep the last (. = [-1] image index) image generated
  1448. -keep[-1]
  1449.  
  1450.  
  1451. #######################
  1452. ## ALBERS EQUAL AREA CONIC
  1453. #@gmic alberts_projection : yaw angle, pitch angle, roll angle, central_longitude angle, reference_parallel angle, standard_parallel1 angle, standard_parallel2 angle, latcutoff angle, longcutoff angle
  1454. #@gmic : Take a map or equirectangular panorama
  1455. #@gmic : and rotate it around z, y and x angles,
  1456. #@gmic : then project it as an Albers equal area conic map.
  1457.  
  1458. #@gui Albers Equal Area Conic : albers_projection
  1459. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as an Albers equal area conic projection.")
  1460. #@gui : Roll = float(0,-360,360)
  1461. #@gui : Pitch = float(0,-360,360)
  1462. #@gui : Yaw = float(0,-360,360)
  1463. #@gui : Central longitude = float(0,-180,180)
  1464. #@gui : Reference parallel = float(0,-90,90)
  1465. #@gui : Standard parallel #1 = float(20,-90,90)
  1466. #@gui : Standard parallel #2 = float(60,-90,90)
  1467. #@gui : Crop latitude to = float(180,0,180)
  1468. #@gui : Crop longitude to = float(360,0,360)
  1469. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  1470.  
  1471.  
  1472. albers_projection :
  1473.  
  1474. # $1-3 Default to 0 degrees
  1475. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=20},${7=60},${8=180},${9=360}
  1476.  
  1477. latcutoff={$8*pi/180}
  1478. longcutoff={$9*pi/180}
  1479.  
  1480. l0=$4
  1481. phi0=$5
  1482. phi1=$6
  1483. phi2=$7
  1484.  
  1485. # Give input alpha-channel to prevent errors with transparency
  1486. to_a[0]
  1487.  
  1488. _fF_zoomaccurate=1
  1489. _fF_zoomfactor=1
  1490.  
  1491. # First rotate the map with the appropriate custom command
  1492. if $1!=0" || "$2!=0" || "$3!=0
  1493. rotate_equirectangular_map $1,$2,$3
  1494. fi
  1495.  
  1496. # width and height of original map
  1497. nw={0,w}
  1498. nh={0,h}
  1499.  
  1500. # Initialize the image using the above precalculated width and height
  1501. # The "" block contains the mathematical 'formula' that determines
  1502. # what each x,y pixel is.
  1503.  
  1504. -input {$nw+1},{$nw+1},1,4,"
  1505.  
  1506. #inside the ""-block it's easier to pre-translate dollar-variables
  1507. pi2=2*pi;
  1508. phi0="$phi0"*pi/180;
  1509. l0="$l0"*pi/180;
  1510.  
  1511. phi1="$phi1"*pi/180;
  1512. phi2="$phi2"*pi/180;
  1513.  
  1514. latcutoff="$latcutoff";
  1515. longcutoff="$longcutoff";
  1516.  
  1517.  
  1518.  
  1519.  
  1520. # Use centered coordinates. 0,0 is in the middle of the
  1521. # projection
  1522. cx=x-w/2+1;
  1523. cy=y-h/2+1;
  1524.  
  1525.  
  1526. # Turn them into radians
  1527. cx = (cx/w *pi2*8/7);
  1528. cy = -(cy/w *pi2*8/7);
  1529.  
  1530.  
  1531. # Formulas and constants:
  1532. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1533.  
  1534. if((phi1==phi2),
  1535. n=sin(phi1),
  1536. n= (sin(phi1)+sin(phi2))/2
  1537. );
  1538.  
  1539. C= cos(phi1)^2 + 2*n*sin(phi1);
  1540.  
  1541. rho0=sqrt((C-2*n*sin(phi0))) / n;
  1542.  
  1543. rho=sqrt(cx^2 + (rho0-cy)^2);
  1544.  
  1545. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1546. # page 114, if using ATAN2, the signs must be reversed.
  1547.  
  1548. if( (n>0),
  1549. (theta=atan2(cx, (rho0-cy))),
  1550. (theta=atan2(-cx, (-rho0+cy)))
  1551. );
  1552.  
  1553. # When the reference parallel is in the southern hemisphere,
  1554. # flip the projection and map
  1555. if( (phi0<0),
  1556. (n=-n)
  1557. );
  1558.  
  1559.  
  1560. if( (rho==0),
  1561. (nY = asin( (C-(rho*n)^2) / (2*n)) ),
  1562. (nY = asin( (C-(rho*n)^2) / (2*n)) )
  1563. );
  1564.  
  1565. nX = l0 + theta/n;
  1566.  
  1567. # For the logic why I'm using inverse formulas, check out the
  1568. # Cylindrical Equal-Area filter
  1569.  
  1570. # EQUIRECTANGULAR COORDINATES
  1571. # Get the width and height of the original equirectangular input
  1572. ow="{0,w}";
  1573. oh="{0,h}";
  1574.  
  1575. # Pretty simply translate radians coordinates straight into
  1576. # equirectangular pixel coordinates.
  1577. X=ow/2 + (ow * nX/pi2);
  1578. Y=oh/2 + (oh * -nY/pi);
  1579.  
  1580.  
  1581. # check whether the pixel's geographical nX, nY coordinates
  1582. # are within the normal range, then further crop to the angles
  1583. # given by the crop values.
  1584.  
  1585. L=nX;
  1586. P=nY;
  1587. if( (L>=l0-longcutoff/2 && L<=l0+longcutoff/2 && P>=phi0-latcutoff/2 && P<=phi0+latcutoff/2),
  1588. # TAKE the value of X,Y FROM the original equirectangular source map
  1589. # and input the value into x,y in the output.
  1590. result=I(#0,X%ow,Y,0,2,0),
  1591. # else output transparent pixel
  1592. result=[0,0,0,0]
  1593. );
  1594.  
  1595. result"
  1596.  
  1597. # If the reference parallel is in the southern hemisphere,
  1598. # we've flipped the projection. Here we need to further rotate
  1599. # the image so it's not upside down.
  1600. if $phi0<0
  1601. rotate[-1] 180
  1602. fi
  1603.  
  1604.  
  1605. # New trick: if you can't calculate the width and height of the image
  1606. # in advance, doing it on a "blank" square image can work as well.
  1607. # In the ""-block of the -input above the last pixel assignment
  1608. # checks if the geographical coordinates we're looking for make
  1609. # any sense, that is, no weird values above or below plus-minus 180 or
  1610. # 90 degrees. If the pixel x,y coordinates would result in something
  1611. # nonsensical, output a transparent pixel. This is necessary, because it's
  1612. # not automatic: you get weird blocks of white *and* transparent
  1613. # lines. Then you can use autocrop from g'mic to do some magic and remove
  1614. # the transparent areas.
  1615. autocrop[-1]
  1616.  
  1617.  
  1618. # Keep the last (. = [-1] image index) image generated
  1619. -keep[-1]
  1620.  
  1621.  
  1622. #######################
  1623. ## ECKERT IV
  1624. #@gmic eckert_iv_projection : yaw angle, pitch angle, roll angle, central_longitude angle, reference_parallel angle, standard_parallel1 angle, standard_parallel2 angle, latcutoff angle, longcutoff angle
  1625. #@gmic : Take a map or equirectangular panorama
  1626. #@gmic : and rotate it around z, y and x angles,
  1627. #@gmic : then project it as an Eckert IV map.
  1628.  
  1629. #@gui Eckert IV : eckert_iv_projection
  1630. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as an Eckert IV projection.")
  1631. #@gui : Roll = float(0,-360,360)
  1632. #@gui : Pitch = float(0,-360,360)
  1633. #@gui : Yaw = float(0,-360,360)
  1634. #@gui : Central longitude = float(0,-180,180)
  1635. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  1636.  
  1637.  
  1638. eckert_iv_projection :
  1639.  
  1640. # $1-3 Default to 0 degrees
  1641. -skip ${1=0},${2=0},${3=0},${4=0}
  1642.  
  1643. l0=$4
  1644.  
  1645. # Give input alpha-channel to prevent errors with transparency
  1646. to_a[0]
  1647.  
  1648. _fF_zoomaccurate=1
  1649. _fF_zoomfactor=1
  1650.  
  1651. # First rotate the map with the appropriate custom command
  1652. if $1!=0" || "$2!=0" || "$3!=0
  1653. rotate_equirectangular_map $1,$2,$3
  1654. fi
  1655.  
  1656. # width and height of original map
  1657. nw={0,w}
  1658. nh={0,h}
  1659.  
  1660. # Initialize the image using the above precalculated width and height
  1661. # The "" block contains the mathematical 'formula' that determines
  1662. # what each x,y pixel is.
  1663.  
  1664. -input {$nw+1},{$nw+1},1,4,"
  1665.  
  1666. #inside the ""-block it's easier to pre-translate dollar-variables
  1667. pi2=2*pi;
  1668.  
  1669. l0="$l0"*pi/180;
  1670.  
  1671.  
  1672. # Use centered coordinates. 0,0 is in the middle of the
  1673. # projection
  1674. cx=x-w/2+1;
  1675. cy=y-h/2+1;
  1676.  
  1677.  
  1678. # Turn them into radians
  1679. cx = (cx/w *pi*sqrt(3));
  1680. cy = -(cy/w *pi*sqrt(3));
  1681.  
  1682.  
  1683. # Formulas and constants:
  1684. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1685.  
  1686. theta=asin(cy*sqrt(4+pi) / (2*sqrt(pi)));
  1687.  
  1688. # When the reference parallel is in the southern hemisphere,
  1689. # flip the projection and map
  1690. nY=asin( (theta + sin(theta)*cos(theta) + 2*sin(theta)) / (2+pi/2));
  1691. nX=l0 + sqrt(pi*(4+pi)) * cx/(2*(1+cos(theta)));
  1692.  
  1693. # For the logic why I'm using inverse formulas, check out the
  1694. # Cylindrical Equal-Area filter
  1695.  
  1696. # EQUIRECTANGULAR COORDINATES
  1697. # Get the width and height of the original equirectangular input
  1698. ow="{0,w}";
  1699. oh="{0,h}";
  1700.  
  1701. # Pretty simply translate radians coordinates straight into
  1702. # equirectangular pixel coordinates.
  1703. X=ow/2 + (ow * nX/pi2);
  1704. Y=oh/2 + (oh * -nY/pi);
  1705.  
  1706.  
  1707. # check whether the pixel's geographical nX, nY coordinates
  1708. # are within the normal range, then further crop to the angles
  1709. # given by the crop values.
  1710.  
  1711. L=nX;
  1712. P=nY;
  1713. if( (L>=l0-pi && L<=l0+pi && P>=-pi/2 && P<=pi/2),
  1714. # TAKE the value of X,Y FROM the original equirectangular source map
  1715. # and input the value into x,y in the output.
  1716. result=I(#0,X%ow,Y,0,2,0),
  1717. # else output transparent pixel
  1718. result=[0,0,0,0]
  1719. );
  1720.  
  1721. result"
  1722.  
  1723. # If the reference parallel is in the southern hemisphere,
  1724. # we've flipped the projection. Here we need to further rotate
  1725. # the image so it's not upside down.
  1726. if $phi0<0
  1727. rotate[-1] 180
  1728. fi
  1729.  
  1730.  
  1731. # New trick: if you can't calculate the width and height of the image
  1732. # in advance, doing it on a "blank" square image can work as well.
  1733. # In the ""-block of the -input above the last pixel assignment
  1734. # checks if the geographical coordinates we're looking for make
  1735. # any sense, that is, no weird values above or below plus-minus 180 or
  1736. # 90 degrees. If the pixel x,y coordinates would result in something
  1737. # nonsensical, output a transparent pixel. This is necessary, because it's
  1738. # not automatic: you get weird blocks of white *and* transparent
  1739. # lines. Then you can use autocrop from g'mic to do some magic and remove
  1740. # the transparent areas.
  1741. autocrop[-1]
  1742.  
  1743.  
  1744. # Keep the last (. = [-1] image index) image generated
  1745. -keep[-1]
  1746.  
  1747. #######################
  1748. ## ECKERT VI
  1749. #@gmic eckert_vi_projection : yaw angle, pitch angle, roll angle, central_longitude angle, reference_parallel angle, standard_parallel1 angle, standard_parallel2 angle, latcutoff angle, longcutoff angle
  1750. #@gmic : Take a map or equirectangular panorama
  1751. #@gmic : and rotate it around z, y and x angles,
  1752. #@gmic : then project it as an Eckert VI map.
  1753.  
  1754. #@gui Eckert VI : eckert_vi_projection
  1755. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as an Eckert VI projection.")
  1756. #@gui : Roll = float(0,-360,360)
  1757. #@gui : Pitch = float(0,-360,360)
  1758. #@gui : Yaw = float(0,-360,360)
  1759. #@gui : Central longitude = float(0,-180,180)
  1760. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  1761.  
  1762.  
  1763. eckert_vi_projection :
  1764.  
  1765. # $1-3 Default to 0 degrees
  1766. -skip ${1=0},${2=0},${3=0},${4=0}
  1767.  
  1768. l0=$4
  1769.  
  1770. # Give input alpha-channel to prevent errors with transparency
  1771. to_a[0]
  1772.  
  1773. _fF_zoomaccurate=1
  1774. _fF_zoomfactor=1
  1775.  
  1776. # First rotate the map with the appropriate custom command
  1777. if $1!=0" || "$2!=0" || "$3!=0
  1778. rotate_equirectangular_map $1,$2,$3
  1779. fi
  1780.  
  1781. # width and height of original map
  1782. nw={0,w}
  1783. nh={0,h}
  1784.  
  1785. # Initialize the image using the above precalculated width and height
  1786. # The "" block contains the mathematical 'formula' that determines
  1787. # what each x,y pixel is.
  1788.  
  1789. -input {$nw+1},{$nw+1},1,4,"
  1790.  
  1791. #inside the ""-block it's easier to pre-translate dollar-variables
  1792. pi2=2*pi;
  1793.  
  1794. l0="$l0"*pi/180;
  1795.  
  1796.  
  1797. # Use centered coordinates. 0,0 is in the middle of the
  1798. # projection
  1799. cx=x-w/2+1;
  1800. cy=y-h/2+1;
  1801.  
  1802.  
  1803. # Turn them into radians
  1804. cx = (cx/w *pi*2);
  1805. cy = -(cy/w *pi*2);
  1806.  
  1807.  
  1808. # Formulas and constants:
  1809. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1810.  
  1811. theta=sqrt(2+pi)*cy/2;
  1812.  
  1813. # When the reference parallel is in the southern hemisphere,
  1814. # flip the projection and map
  1815. nY=asin( (theta + sin(theta)) / (1+pi/2));
  1816. nX=l0 + sqrt(2+pi) * cx/(1+cos(theta));
  1817.  
  1818. # For the logic why I'm using inverse formulas, check out the
  1819. # Cylindrical Equal-Area filter
  1820.  
  1821. # EQUIRECTANGULAR COORDINATES
  1822. # Get the width and height of the original equirectangular input
  1823. ow="{0,w}";
  1824. oh="{0,h}";
  1825.  
  1826. # Pretty simply translate radians coordinates straight into
  1827. # equirectangular pixel coordinates.
  1828. X=ow/2 + (ow * nX/pi2);
  1829. Y=oh/2 + (oh * -nY/pi);
  1830.  
  1831.  
  1832. # check whether the pixel's geographical nX, nY coordinates
  1833. # are within the normal range, then further crop to the angles
  1834. # given by the crop values.
  1835.  
  1836. L=nX;
  1837. P=nY;
  1838. if( (L>=l0-pi && L<=l0+pi && P>=-pi/2 && P<=pi/2),
  1839. # TAKE the value of X,Y FROM the original equirectangular source map
  1840. # and input the value into x,y in the output.
  1841. result=I(#0,X%ow,Y,0,2,0),
  1842. # else output transparent pixel
  1843. result=[0,0,0,0]
  1844. );
  1845.  
  1846. result"
  1847.  
  1848. # If the reference parallel is in the southern hemisphere,
  1849. # we've flipped the projection. Here we need to further rotate
  1850. # the image so it's not upside down.
  1851. if $phi0<0
  1852. rotate[-1] 180
  1853. fi
  1854.  
  1855.  
  1856. # New trick: if you can't calculate the width and height of the image
  1857. # in advance, doing it on a "blank" square image can work as well.
  1858. # In the ""-block of the -input above the last pixel assignment
  1859. # checks if the geographical coordinates we're looking for make
  1860. # any sense, that is, no weird values above or below plus-minus 180 or
  1861. # 90 degrees. If the pixel x,y coordinates would result in something
  1862. # nonsensical, output a transparent pixel. This is necessary, because it's
  1863. # not automatic: you get weird blocks of white *and* transparent
  1864. # lines. Then you can use autocrop from g'mic to do some magic and remove
  1865. # the transparent areas.
  1866. autocrop[-1]
  1867.  
  1868.  
  1869. # Keep the last (. = [-1] image index) image generated
  1870. -keep[-1]
  1871.  
Tags: gimp Projection cartography g'mic map projection
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!