From 36708be6f77518e76220c5af0848708bf5762e2b Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Thu, 25 Jul 2019 13:51:54 -0700 Subject: [PATCH 01/30] add a new module for manupulating URAT files to the API --- doc/api.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/doc/api.rst b/doc/api.rst index 7ba2d14b3..99ae31bea 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -88,3 +88,6 @@ desitarget API .. automodule:: desitarget.train.train_mva_decals :members: + +.. automodule:: desitarget.uratmatch + :members: From c9cb5901b997b2d2c536c0fba72109639f3d9dd6 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Thu, 25 Jul 2019 13:52:51 -0700 Subject: [PATCH 02/30] modify the URAT-provided FORTRAN code for parsing files to work at NERSC --- py/desitarget/fortran/ADM-v1dump.f | 97 ++++++ py/desitarget/fortran/v1dump | Bin 0 -> 31064 bytes py/desitarget/fortran/v1sub.f | 541 +++++++++++++++++++++++++++++ 3 files changed, 638 insertions(+) create mode 100644 py/desitarget/fortran/ADM-v1dump.f create mode 100755 py/desitarget/fortran/v1dump create mode 100644 py/desitarget/fortran/v1sub.f diff --git a/py/desitarget/fortran/ADM-v1dump.f b/py/desitarget/fortran/ADM-v1dump.f new file mode 100644 index 000000000..751ec28a7 --- /dev/null +++ b/py/desitarget/fortran/ADM-v1dump.f @@ -0,0 +1,97 @@ + PROGRAM v1dump +C +C gfortran -o v1dump ADM-v1dump.f v1sub.f +C +C - loop over range of zones +C - read individual URAT zone file (binary data) +C - output all data to ASCII file, no header lines +C - user select column separator character, no empty columns +C - for explanations of the 45 columns see "readme" file +C +C 141120 NZ adopted from UCAC4 code +C 150119 NZ change to "v", bug fix "238" +C 150126 NZ add URAT star ID zzznnnnnn to output +C 150129 NZ change format item #6 I*3 -> I*4 due to neg.flag (no GSC match) +C 100725 ADM Adapted for NERSC and hardcoded to skip interactive steps + IMPLICIT NONE + INTEGER dimc, dims + PARAMETER (dimc = 45 ! number of URAT data columns + . ,dims = 45) ! number of separators + + INTEGER dv(dimc), mi(dimc), ma(dimc), is(dims) + CHARACTER*1 csc + + INTEGER i,j,k, jo, uni,uno,zn1,zn2,zn, nst,nsz, idn + CHARACTER*45 line*223, a1*1, bfc*1 + CHARACTER*100 uratdir, pathi, patho, fnin, fnout + LOGICAL eozf, bf + + DATA is /11, 21, 25, 29, 32, 36, 42, 48, 52, 55, 57, 61, 65, 69 + . , 73, 79, 85, 89, 92, 95,106,112,118,124,129,134,139,141 + . ,143,145,147,149,151,157,163,169,175,181,186,191,196,201 + . ,206,210,214/ + +* defaults +C ADM the input/output directories. + CALL get_environment_variable("URAT_DIR", uratdir) + pathi = TRIM(TRIM(uratdir) // "/binary/") + patho = TRIM(TRIM(uratdir) // "/csv/") +C ADM the file RANGE to convert. + zn1 = 326 + zn2 = 327 +C ADM byte-flip or not (big or little endian). + bfc = 'N' +C ADM the delimiter for the output files. + csc = ',' + + IF (bfc.EQ.'Y'.OR.bfc.EQ.'y') THEN + bf = .TRUE. + ELSE + bf = .FALSE. + ENDIF + +* prepare + jo = INDEX (patho,' ') - 2 + uni= 11 + uno= 20 + nst= 0 + +* loop zone files + DO zn = zn1,zn2 + WRITE (fnout,'(a,a,i3.3,a)') patho(1:jo),'/z',zn,'.csv' + CALL open_zfile (pathi,uni,zn,fnin) + OPEN (uno,FILE=fnout) + eozf = .FALSE. + + WRITE (*,'(/a,a)') 'begin read file = ',fnin + WRITE (*,'( a,a)') '... output to = ',fnout + + DO nsz = 1,999000 + CALL getistar (uni,nsz,bf,eozf,dv,dimc) + IF (eozf) GOTO 91 + + idn = zn * 1000000 + nsz ! official star ID number + + WRITE (line,'(2i10,2i4,i3,i4,i6,i6,i4,i3,i2,4i4,2i6,i4 + . ,2i3,i11,3i6,3i5,6i2,5i6,5i5,2i4,i10.9)') + . (dv(j),j=1,dimc), idn + + IF (csc.NE.' ') THEN + DO j=1,dims + k = is(j) + line(k:k) = csc + ENDDO + ENDIF + + WRITE (uno,'(a)') line + ENDDO ! loop stars on individ. zone files + + 91 CLOSE (uni) + CLOSE (uno) + nsz = nsz - 1 + nst = nst + nsz + WRITE (*,'(a,i7,i10)') 'numb.stars/zone, total = ',nsz,nst + ENDDO ! loop all zones + + END ! main + diff --git a/py/desitarget/fortran/v1dump b/py/desitarget/fortran/v1dump new file mode 100755 index 0000000000000000000000000000000000000000..025d0da9d2874b0a0cc8b9d6015f7f93d836383e GIT binary patch literal 31064 zcmeHw3v^V~z4w_+0zn`XKwH$Rr#@P+$s|BPf}oSY;KU{xktkMh7?L5GHkrg^2Ei+e z#_44`j;2*xRzDwHYrDLD)>E-(m{7kLCqo(*HLU4RO=(3oU+ zjL#Y87!!cW$Keqa#xkYjQKeePg^JIKm~spABw}HnMQ$Kkxx}<^Z2J71%rU5t|hptPPBF;qo zDIBNaD8lh+98+=3z`?T8os)2|T~NnI3ED52Hk<%q)-CI7ic&OL#nV-M5@KDJlaW3J z$B8&j#lgTb>O&wt;^`*v={)Qrx1)gB@K@&G|7;HYA7E5AJ!Luc-;slVN)CR14*Wwo z?0G7O{*5{Cr{%!kkR#vqIq-{e=%H^+9GPKV4*s?r`X9)_AA;g+{(nCQK9WO!T@L-n za^Uyo(ErUG_zQF3U(JCx(C?eqw|7Q517Y^|3d8SjZfg(wyCQ+kh~IDcHPug+uJyrA ze`hez>c6ft6bag>jqTyaKm=l)p>VSe-xVQ$QyQkdBPXA~sXc603-R`jVAv+!*xuF= z=nRh6VA&T6M}p0;Gt)%q5~$2PEWyT>Kxd#al0|epi_IeFZES7t$}l?=zBU@_OhX1) zm`JES+?6IxkNHF4rr^dj^f>!1BUAxx+|5uO3~vZ^wujq-;fQ}jpfeO$*BbPPP@z^0 zc7ZF<*_$SUM5AeJ^mh?ChtRdrKvP>F(n3Tz z5si|+E!c)k4MdG?9Y$NAwY9yGP+R*3c*ajg)zaDmp|G*OGZ-{3s;#MB>YwYL>z;31 zwEWT)D_2xq;y(}VHxIK;p7_fbXAKe0w4EbGqXfrvhj9)P9Q#ci(NjX0cJhroNkkyM zVI19=Mn@MF^5k#+>o+hB~5AmyPc+S&0thM1e@9NNE!|S=1v<@3y z&&|YdwBflf&|!-W&vkUe3 z5IXF#;knMyq1T4znn#Cz8=h+;9R_T8u1jgzsGxL#2oDZ(uz4j&#vIOk!) zY}&&E0)BvS9^qaAf15Cy^6+B<{sv(-;oi@)mYer=?FWldkNaiPcX_3iDsfKlJ$Gd;%8uOayEde<4t z45MM!;ALwK-x6_-|1t9ET~x!SGj+|=M2;O6d~5bg{yoGs_-_8Q67cQ)!0C&R`1TH6 z=yM$OJ^w-EBum5OlmkJckCvk;`I51bzm1Ktt$(?x>Pa@9Dv|4xl-bRbFMjZ86D1eup`Cvx z_km(x@@^rR+%Z7;Wb9c+@#l9p36Z{`HSkPAV7q$@ecN`kAHW$t*Rp-Sc;C=(p{FK( zz?Zn2@{%!;q>|*O!m1}{A;YR&GzU84`)cE_4=uJ7#(Hmb`4TaaF;v}oTq{V{7y6Pr z4oibxK~x+6Ju`bAnI-Ne!q>Q8TG&0>jzLS-9>rHm`GAda#KOoljXu)94ah#O|=v%b4124r^!-~tS^AOdTB{A_ESX4T@SD% z_Yz?hKIJaOz7{;nU6(VN(_Mnm=B@^j&l^L(d7tHACdxjc81@|iq0c+C91W(=``Sf_ z4tjri$wBX-64y@eFBe}GDHQRh=%haH3v5+KCL)2DjhVTr20R}=k*y95n5}Z?1Z37d zlo$6N?u)GkBZs~|?<*W}`eKXWFC>hFcPO#uwRr3(m=oe$x!-k9uWM&-tk)I0hJB_l z*8kVBu?5iI7aKro3B-^$2N7R9X43xHESyDG=!+RhQgdJI0GJt|p%BQWLPIf99CZ6) zy&xO)LKl;4iW5N;9j|-yVZ&A9JsjVMCMbNG6#S|Vo(3OGpxEv_DBE)>bW_3i1R=@N zXnH>AO03@|M)aP`lvsO;KrI7O3x{QNe5q^)T?n1*4^>aHtW~>2iQ|28sHn>3D_(`A zd^G}l`(j6Nret62A>@OE_63%n3GXX$;WuSHtvM51eX%KUBzebGN3}%K!wr^)Ic+Kr zV}_DdL7YwP{~AsCh_lq5tCcg>u`-fLS$gSXSx_-rEE7~!8GS?FC79VlmZnT<`#2^Q zW}4KcO!|c&Bt?@`)owHChbb}oYdTEoi!DS0XLm`6)RI*fxD40*x?Pu`(fE>2ah^v> z=~`d1;RTEil!A{~*Huq$Mea56QRtL4{2&$+2~iR38$-}8drWA+Q|(e2IR@gf61JO6 zpKzNInmx?G-fMm*eU-D<^jHq$&O_E-^ZYRs*n3o#mVKi8u(SKH%pFJ4!lA-itYTb zmQS21pP-#*>pf5WCyh)BzpEk$m&;;qfa`2U_%%S;V7hx9anW*ye*&AvTXR(^W26$I z#feL}L&}Q8-5jtbl1+@zIL&LI8Ay8QxaA=(fd9cr@bz~s zC+z0>anA&gM~tScr0-Ox5TmWFiCizK!j$oiW3`}K-vrG)%hNjxn<4Xod zuC7|WWgt5enq2k69Q}C_mMy}f*AWvUx8+bhMbYE;cQfEu#AXN$EVWi~j z(t{MYLHc5MDV#4EBa(x`lD82{aQ9izv8Nu=B;=zeu?!XEP;^_H`W&e#Ef=xcA#FU( z)6gy|P3)jYWO3Enrwh8p+6Vei1P_Ljxc8NDWtj~fC^UVCA$UY+p^xIRQKGQALBxz? zi7+b=$)ZQBqRZ~YQmUm3lPp*5Th`L$6)s&~iSI)(X{uEoF@NUU9K#WN&yC74N1!0L zW2UCWn46nmW?g;Z!^@Y9-M|vXW4mBfpQx!V9$i?7a(dFGxe>-B1V5&_TVo_<=H@2S z@@6k|$`%)o{gj24j=M=X?gb!Zn_<>i^1;I!F-anU6!Y9tU^t+RjU3%eG$@-jCaf8r ziAFlkh&ePuW)-vQi0w(L&DlKf`wL7|&5xJZd*Z4Vf+H4zem+Prv*dGAnci<@3hVw+ zrn8Z0;%?C^6QT!OMHDUL&6G{d`it4Jp7&2?{pg>qV(k{?eL$5rnx5&41T)K9@lR%2 zY0H%RCVO_NOpWwR1KBd2|4}m)mw~C4d+VQA(z-3WcG=q0nr2eU9y)lD;Ph^J(i%Q^ z?N9!pAAp3NjH?s(a)iazUTt#iW3iVb7xM~HXzU)X7Z=pbDG?LhT&^f8n2ESNtlw~| zCUy7L*T!GnV+vdPhCIw!O3BNYKINEymHB=)jrfcrepwJ*vG*~*+Xlt&P}$7~I9p&E zty-fdtDBEL1`OtpjVwWYM~8@H>x(DX0;-M2R*5KCPg)|zC}CiQ824iK^^VB>5N=}l z_St6<$=5wv;=1EjmZ(pRD00r=IIvD4N0F5M!$4n`{o^s^+UE2rPSjmJmlvdwg_ z91$egB@1jL?Ukf^GznSgCj31~+Mb1!kfa@&L^r60-yum)WFa+360Uuu8MeF@Nm6eX z(ixI;UlvlnB<<8Bn;FB>8{f%7dQOsdYm&{3A4$@#ETnCcv`3R{d2N=Yr?QaNNz%ib zWXtO!NqRI3sYH?<$wHbYN!Y58KD3!}Y*^GXHqyqCUY4YOO|mtBUP(HrNj3|=D@lj4 zkiI5K2eOdbC224VX{98+l!a6&NiSw0oi0fOnq(`oQ<9#~LK+$p-uOip($6L7@hqgr zCFw_5NZ*pAA7&xlK$2Vq#Mu-2#Bx^HiH#R@Q!#JIsZL;%STV-M3Ix0y&>65Lz%_}| zFL#M$ZS)O&`HhdA!|NX}heGC%tow16;@vArPiG#29$%0=UDM^gc*ov06`hG3NCRjfy3x_a3It5$kmCcBW63=@*!u zD$~78J7oGkrvFpQ?q+&ergt&@8<~EX>E~toQKtWoOh3Z(lQP}UbVQ~PGW~BdeTeD1 zW%>ZqahV=udb3Qw#Pqc?{UXyrr2D+$*|qM?qXpJnu^1PQdW!_@!Ewmv7$E(-k2P{M(5ZmpDDY$gp;o6hpCIYjo@MJj-J%8pnF zk%#Rfi-CC5j_?5Sh#j#EI`-Lh_<-1NN7N=7`YmtOL;j#$ekBlx?1)uB9Izu+12Jev ztOepFJE94Q7ww1^*f3zXp%sYd?T8K_eql#MfOy=F*a*ar?1&qH_@NEavjt-ja(=~@ zbI+|3_gg!TD-`e@w&QM-xL57C9*KL+j=M|ZUbo}6N!*YfcaOxqX~%7sxTAL5gA(^= zJMJOkc6xtoFWVy$_p}}N7;rsugQn^Vy|ZCr2*t{DkRxk+$AF0L?q`JCT!_>n6urSa zU@b!9zx26k2Ejo_#ep$l(C7Mrx8LFtOV!>gF;dsWe_Qk3i%RFRgxTjEtZDas(nh(Rn%yP@OA zQ$5|t>BzUZYd8dB<;r*jsX+?iZVA@Ia%IdNZMiZ&46<;@4^U~w@}b&a*J=J9d=Djq$I4Kj?4bpy3=Z+~sP zt`}Yw+h5r7Dp-&E&tkKz64=As6eG4?W6uJ1W*>`=vU_{w41!d##0~QyJAgyGi>o}w zo|!n&s%yjo-FrA4>yuc$(Hak$*snJTsmX-U;b3FCf@(vn|t4%9`D~f zddBR1@quJLHYU8ig0r8TV}l>cSu8k<6U=%rVN~wBaf(P0mb6|Q4OfWdvS*_Fd*oya z2-e<&k+NiFyjbmbb|3xFD&1|e4rZx3c!_INxGo)gKw>Q)UWdmMgb%~Ng%5GrWNQRh zYuXAy8^u-C`Pc`&;~TOupb4<|ByC|CD+@~ICN^dgBA zlY`h=Rdv=dDiVV!;($sHCnYm1;!^D&S|qYCwqrerRKLpiBbZsApJoxdK+}bFMA=MA zz|8u5H2~rvgtcq8Gi^u^ zralYg+pg7nEi`+D7$ZT&miAwIMW?cfiOz7h5L=a9Oyt%Dj_c3}$LVr?m6HjEF^+w5 zY~#SRm3_Z2HWMRSIT~GG%tQ)0km8sw$2RT+@&7Hx$axCGk0G>*EEm1AV+V;yu8`=;G)Zecl6b74B;F3Or5j9-SG* z{r_I#Q8jEp&{hcT1xK*OM(u1V?;)L|n(qB4tEu%R%ip^(kQ#fgRMj+(6_H4awmiO? z3Q}UMrVU^o9%Co4j7RkpPdyz5>uMT3=>6RnQuI_CJ_UqQaMz8pyr+utUMc*|n?0Nn zSW(X=ywB>I1u9juhbhWa!9#eei^H50U=2$AHigG9plf!sL|HXEz?^uR+9)WpX4{e0 zHG8K>qh`<5H9G)(q6Yh!K+R&8s-HCoO9lSfah+)wYbhc2)o{tkRRZp4X(BJn3Ad9|UUkC>6X67`0>nl5LC?hUp45b} z%SoG_8pZukEQLozboY>m#G_i+q#XoQ8&3`pDIWtAujrLZq{2)nq?*TfI0dx0RncL)z037^QP!V;6>b(ry!%#Z02%$)JMXNTXw z%-EBN9#7PbINo~4buPv-jH`H>=pA$+roS)JtB@|~kJsTLwPRE2E>X`}ewb zJ`;N;x+qaM>T7%khQYs!Pe_2f@tN9WzDwbgn>d!x8`X&^H^*;#Y6MSWh^y--dxSYh zVMV+S_hFOtRbrB*tOCl86LLJKG$bm1n8wyU>`2^}q{7PQT(|#>Z5gD{FkKa*#(@h~VY$p#lp$Vmvo6@RW0_L47E*~PcF-mzT%8gTK4PxQM@1=B znpyZrh|@>Ii!l)NfR@EIBOl#0$q~^dWhH)O*Qg9K3rX8QzH6xLzs)tDv>LRU?lwGv zleOgxedLyN({5<^=v|+RKdkE;sjkmwg0fq8^dnn$!;^`+qvzKRx$b!2!@GX;rc@*L zCD{R0CMmS+_AJ21Z`)KR+cpb`Qt6RGbm0iQ@NgLOd!lYA{@fSDu=|C&qm|FySe~eR zjipU;C^;U#ZO4G%O5uip<4PO`Vy8s} zml4^VU-xWuiik&Cw;jU+J4ko$b;RrX6VZN-#nE>w>)!Ra?my^p?K~Dc=8F9Sqemk8 zZaj*)WQS-C7o>*Q6l3)KL=;oOF-!#`aOYfz0IEw)DVO-#9cq=VyS;LEO7Zcfo=xc1v|3$m+@Nv~AJ)G0- z#md0sjbph@P5e~MvdZG>xU%yJ>|(v^7(N}13>LrU=-&LAabu}!Q*th414K(p!?6LP zH4Xm|F_t-a`2SXPH#{wSWmSBv>WZ59v8sle_^+x~)+Wy}T_%1KkSt|aXFuk5{-#C!gvB8?1B5B!)UcvzrGHt?in4W~aG-NjSy^M(hBC8iSzYOdIZe^F z4)=OvR-hzM5}N0pR}z?AcD)-6v&w*)y)M`s3Y+{fqq#oR8Z?)fvrOXMZnxPUjdVmK zW~ALTfk#fID$7}OLvzYY=7uVuv4rOZ;$MPuODcexOI!s6vj@;2IcHAEJn+p6%`aIH znp-j->G?NJ^0uoF5>%VD#xZ^ADg z8zb=7?8ZRY47W$%s&JDTX$cCSErOLDt-(N7(8RA4Asr50XLjIsgNtSbH-;kSb)iU$ zdFIA5XKyrDnZlxAXJ>n-xfs8146!^HKqA;=Hig<6%?p;874yZoR98TRPhWv9jE33O1XIW}&`n@H52{i0gRydG3nYv#RD;2#c`5J&zBf ziJeR{a;Re^rGKie=F%1p786l7CTOLo_`kr?$DI~Y;5c>V&mDdu|fxa+lcrvfNvl+ z0e|B~>`Wp4IpP+?Zd~keyQdlPLx{hLxF7M$zZo0j;xUYSuf@2?e-pQY7b1QN&jhSO z+>Racjfl74&fqr09hd-jAwG;4Kim|7dkYMX>sJ_#jm3^rr%WvDaZD^G{xdjM{ok>% zbIG7Wc7DS+8+7vfo*5e(fK5hG@$#aX7rQ23SGdKv@Rak;pMThJNXo&=5yYD0vV{M7W6A{yZ8pg7X6ajoyq*$3c4paLk=*~o(Sj< z<1X_r?R?JXNEYS23-~pt`^VGh-?!*b0RIzpDCrGEa-I7fW!tRN6MV{88g;ho6 zRYkL^i_Gexnc%N3D&%&gu^4S{5PA@O5cEAOl{O0rMRv3{zxeluz-96z?!Gv)$4?ma`t}qF z8UuA}MT+Y;22IbGps_K}BU1nQ!B{&J9;bwHb3ugi2@>X#kYTN&bE(Sku(FrGVrG~r zBLk17qy=>ufft#q6?i3PY*z_=y};j*G3a-4^i1*ySmQS;SeHZR$4dbQ4Py&LlxL>m z$7Pv-wH;W-3izc$kC6NCI_dY=W#+eD)$;}we^tfbRPpyzyi>&oRs5ogUsv%l6;GHd z%k?Q0e@?{fpSEzWsiZ`hEt1AAciod7gohm-4;ulr?x{8mfc!KH|pHgvVd2_*R zM{GN6->2a)aqt?5*OuitDsb>!L%b$if};wD7e@__OK~*dxC+N=9DW>M#?gUe1CAST z+>GO^I1)ILIPSplbsTr%xEsf}aNLXIJ2?Iw#}9Dqz`^@@Pvba%;~5;!;W&ijMI09` zUAoAeb$P?(UbB`j=bLR-xEJ!}NCJ;H^368SbI)^^n{&(OEG(ZhcfL7mDPG%T)`T0~ zvyH0SmCO8)`-;CNyuKakM7Sc@*@YKTEt=~tUu4t;8e2l);G)L%&d?%bWdLuZYFbn= z_q+vx;JS*^`Q;1emsZSQxUjTxVbh$_IgRVe*98J|&zo1VevxriFwnUO`BpB{n$z?1 zd-1v`TPQ(Vc1hX8^m{4k|7?NdN}l1yS@+0vuEm1>h(*`q8S9uqk7xN39Cw!X(A!c! zUwgvPDbC|hGQZEE2$J8c`hyr1nY9Vr~6J0LA+p^NbZp6*{_!K)4g!Df?DZfx!wcbWB6B96wj&OgFv( z7L~p85%A=_z#S=NY`!$FfbSG^EY+ec6occQKj1upk489_E*B_(e`RhlB!KiqC);k7dGZAS(EwBuZ8pEca273^eE#H`85f?MHxqB{{G3rULr6no!-VYmkkiCW zou9H4L2D|50=MJBlVuj4Kbdl#ncu2D*?}x36j%~@uy10@5|k4)@o6$m0=7`vN2U|v zw>VHP!JoMcx+1e}}co>jo+iAlrZ_*>Cb!+}6X{JY&{U7Na^f*oCD zc#*1q0UC^WpQ8J6l2OP!-S3?jIGn}q>a#AWJKI@gMmKD{IuGvyahxS|pGI|=x{Xjc z6md6Fjj7MLIJPV%GIkMmv0bYplegOu?6pjY>#%8M#I

IN@S0_G^Ojo z+xV5z#@0YrOKC@EdjvOru%ZAVGlxiPS7}>&Q?PXjt3MQIEe&oAHmX9i`2Lzp@fOIY zP}q;A(_~n#;rlqHw}g{I*9Rk;I)cXdDp4)F3mpE4ffr$R7oBoWlCPjNC=L;FDtr)&7xk7(>cG9Mkf_W5Pb#LqW+|#GNAZD%pISju}>bX4Mc*oIjzn0C6wE~B+QYY70>aWpByRMZVLDywk*IG;08lAN3 zx@?@Tc}ifRiZR&VuhUmp%2K{isI*AkuN%`yy*|}rj~;IIxTE=W%=>f9RK#G4 zvtfKz#d`g$ZPxZFvK3|8qwTm@%?8?zscKzg34Vlx6JzG&x| zz4B)|eyU;G*7=&|#GBFedO+hCI}Mk0{gRF^HuOq89Y4)_zkE7=y0NR;V@HX!6K`SH z`#W}`6IUF1Und=Z65gA?){YWsC*JU`-+yZ-I=RbIpN$`_`Lb;IlZ_X$;<@e}ADq@( z+tX=IaW^s@N60Cz__N_pACDCLXXL<}Iq*H&O^-)lKBMpl5NFf>Yz}-bfBqW$4!qI5 z!s|gQJ!I4WN5L{b^NCW1gm)F0jbj!e*1ih7tam=&wcT1_HaKu8$P7;+-}3w z&mCn_k9bxCdhQ3_M4TBO1%4W~5ehF`n@0Kx@yJ2XQ~8SeGoe3?RRBi?zR1v*fqG?j zJO@9&>ypiHZ$Voz#y7rh*dBV6em&3X_0oAlPm%GbnI6F&P(sUqrye~{YF+g?^6kul z|2L(d?@?wD&kG^uZ`RnJFH!BeQcw-!`#JRN%Ypxc(({{AkI*2Vy+P_N;8`!3e)xM1 zJ!ceT_h+rZ7a9Mm?Brep!{s^nj|o4IUl(1Qga7UvcmtPoU6ClB>})ief}MEuvWwqk zh_vxj3*jI>VSukV_?uhX*Wp_fO_BD_E`K1p(cs4+T7!{blY2qsy!qLg#ATdyffwoA zgwGBH+Jb&O{@Jz(A~q=fT9B-1Pm0!6)m(y4WQ4lh^X8YIXQ*#GxaU*?K5w4&EM?kN z=kgU*bzZ;sl4bn1M!J4nYc9L$lB&9zrRfy0;Q$o21pHo~Uuy6zTVePws=c(js@8w$ z^5vI%SNd00Ro8l{ReV;#x?r_^|3Q9>Lwx_jq>okBuz~!uMhbYM-B=3%@ni+d;uf}i{)l+ ziHJ{>q%&rINdr_{Wy;TQq~r7xwHaO%A2P|%P4}xUl0etombNaovGk9`_?x=g{VjoT zQ!92*(sy%Go?(rvK9Z4x)nLUgjr>eVy16MWmMgF)#Lv8>sS?j8TW!t1uB%HKo|5m1 zSf8Is6HRegIatn?pJGX4m)|Qvqp1+Bh95M^Wih@Fqdo?crhyHmGmtWNW!+LeE?6y2 zH4T2AB~3PqvC17;Ol8xan^Di~O*u5CeaR;WPuiz+(s&Ga*QPdXE+UR}%D6?RXkurF z;l^{}L3dSk4K^X04L3F`-7VNeG~7*_!Vs2mq*GEh$QzYufFHEZU~7N^D%H^%G2Ei{ zxRG=>w*!m>H{vWd?cAO1Vi0i$Th#d8(gY}>kPLU&v7$N|T&% z|7+({A&#XruEAtjnCB6Rg$kjsuOwCS>NMlCaPS!kEwA@;^!M;g(Ch)rjI=Nx5T<%d zeq){o5WfGS6dn*z>ico>?kcbKWSD~pFaB-U`FormLGWO`!JzGD9t>JuLUP^3d;eN~ zt&(4>mAo#e&R@%G{IyEHS}D}~ zVLeKIyHZSf?D|;wYx#cxj(%0{vi8l&rvQnNxx=IDT~g)q6?S>OzgDm02S~JqVwKkM z&4BIldLPc@1`2{MuY^*t#(xb2yS%=RC|B}_%C)47b-o(glOeD7ee*Ar9LkNF-D>w5 z@l7SK>rY=-6)X8-)s8iPX8!kO$eYz3!D^}sy-Pr;Q>OprK256q?5g$%=3PpDi#AND z(0)^(u@ecpN}H+o9rbrAGx;^F<9(2I+2r;9KpQebIU;F{%YstmSlk7eBY>>)IETd}e(!FMBuwfVgry5@{`eK?YvS>-b&Jb=uB4eT9U& z?Z_;=?`OeLWV@2jRIcldJPebNdA9twtJ^oejP{$!dr}Vhp$4hHJ3}A?G|tK)|CKAH z{OcJ48K6<7Ytyo0VvRXo?W)L$#gZpxsI7;^L20X%m zjI48}|0$c{V!^xi8Yv>a*oEZTI2Pg1?Ox9x93L}^2n?Q2Z4r_2q@K7WIL_8-z-N=c OIV=S(&nR35X#77V`jg=R literal 0 HcmV?d00001 diff --git a/py/desitarget/fortran/v1sub.f b/py/desitarget/fortran/v1sub.f new file mode 100644 index 000000000..6a5dbffdb --- /dev/null +++ b/py/desitarget/fortran/v1sub.f @@ -0,0 +1,541 @@ +* subroutines used for URAT1 release utility programs +* +* open_zfile : open direct access, unformatted zone file +* getistar : get all data items of individual star (single record) +* nx_byte_flip : apply byte flip on index array +* flip2 : flip 2 byte integer +* flip4 : flip 4 byte integer +* valid_range : restrict R*8 data item to given min,max +* get_zone_range: declination range --> required zone numbers +* get_ra_range : RA range --> required index for RA bins (900 arcsec) +* as2hms : convert arcsec (RA,Dec) into hms format +* sorti : sort 2-dim integer array +* +* 141120 NZ taken from UCAC4 release utility code, adopt for URAT +* +************************************************************************ + + SUBROUTINE open_zfile (pathz,un,zn,fnzone) +C +C input : pathz = path name for zone files +C un = Fortran file unit number +C zn = zone number = 326 to 900 +C output: fnzone= name of zone file opened here +C +C - construct file name from pathz and zn +C - check existence of file +C - open file if available (direct access, binary) + + IMPLICIT NONE + CHARACTER*(*) pathz, fnzone + INTEGER un,zn, jp + LOGICAL ifex + CHARACTER*80 answer + + IF (zn.LT.326.OR.zn.GT.900) THEN + WRITE (*,'(a,i5)') ' invalid zone number = ',zn + STOP + ENDIF + + 51 jp = INDEX (pathz,' ') - 1 + WRITE (fnzone,'(a,a,i3.3)') pathz(1:jp),'z',zn + + INQUIRE (FILE=fnzone,EXIST=ifex) + + IF (ifex) THEN + OPEN (un,FILE=fnzone,ACCESS='direct',RECL=80) + + ELSE + WRITE (*,'(/a)') 'can not find the file:' + WRITE (*,'(a)') fnzone + WRITE (*,'(a)') 'please enter new path:' + WRITE (*,'(a)') '(exit with "x")' + READ (*,'(a)') answer + IF (answer.NE.' ') pathz = answer + IF (pathz(1:2).EQ.'x '.OR.pathz(1:2).EQ.'X ') THEN + STOP + ELSE + GOTO 51 + ENDIF + ENDIF + + END ! subr. + +************************************************************************ + + SUBROUTINE getistar (un,nsz,bf,eozf,dv,dimc) +C +C read data of individual star from URAT binary zone file +C put all items into data vector I*4 +C +C I: un = Fortran unit number of zone file (assumed to be open) +C nsz = running star number along zone for which data are requested +C bf = byte flip requested = true or false +C dimc = dimension parameter, max.numb. items (columns) on data vector +C I/O eozf = .TRUE. if end of zone file is encountered +C O: dv(dimv) = data items of star (or zero if eozf) +C +C 141120 NZ adopt for URAT + + IMPLICIT NONE +* variables from argument list + INTEGER un,nsz,dimc, dv(dimc) + LOGICAL bf, eozf + +* data items of zone file + INTEGER*4 ra,spd, id2 + INTEGER*2 sigs,sigm, epi, mmag,sigp, nit2,niu2 + . ,pmr,pmd,pme, jmag,ejmag, hmag,ehmag, kmag,ekmag + INTEGER*2 apmag(5), aperr(5) + INTEGER*1 nst1,nsu1, nsm1,ref1, ngt1,ngu1, mfm,mfa + . ,iccj,icch,icck, phqj,phqh,phqk, no1,mo1 + +* others + INTEGER i,j + +* check + IF (dimc.LT.45) THEN ! larger is allowed for extra col. later + WRITE (*,'(a)') 'error need dimc >= 45' + STOP + ENDIF + +* initialize + DO i=1,dimc + dv(i) = 0 + ENDDO + +* exit condition + IF (eozf) RETURN + +* actual read of data from file + READ (un,REC=nsz,ERR=91) ! 80 byte + . ra,spd, sigs,sigm, nst1,nsu1, epi ! 8 + 4 + 4 + . ,mmag,sigp, nsm1,ref1, nit2,niu2,ngt1,ngu1 ! 4 + 2 + 6 + . ,pmr,pmd,pme, mfm,mfa, id2 ! 6 + 2 + 4 + . ,jmag,hmag,kmag, ejmag,ehmag,ekmag ! 6 + 6 + . ,iccj,icch,icck, phqj,phqh,phqk ! 3 + 3 + . ,apmag, aperr, no1,mo1 ! 10 + 10 + 2 + +* byte flip on I*2 and I*4 items, if needed + IF (bf) THEN + CALL flip4 (ra) + CALL flip4 (spd) + CALL flip2 (sigs) + CALL flip2 (sigm) + CALL flip2 (epi) + CALL flip2 (mmag) + CALL flip2 (sigp) + CALL flip2 (nit2) + CALL flip2 (niu2) + CALL flip2 (pmr) + CALL flip2 (pmd) + CALL flip2 (pme) + CALL flip4 (id2) + CALL flip2 (jmag) + CALL flip2 (hmag) + CALL flip2 (kmag) + CALL flip2 (ejmag) + CALL flip2 (ehmag) + CALL flip2 (ekmag) + DO j=1,5 + CALL flip2 (apmag(j)) + CALL flip2 (aperr(j)) + ENDDO + ENDIF + +* put items into data vector + dv( 1) = ra + dv( 2) = spd + dv( 3) = sigs + dv( 4) = sigm + dv( 5) = nst1 + dv( 6) = nsu1 + dv( 7) = epi + dv( 8) = mmag + dv( 9) = sigp + dv(10) = nsm1 + dv(11) = ref1 + dv(12) = nit2 + dv(13) = niu2 + dv(14) = ngt1 + dv(15) = ngu1 + dv(16) = pmr + dv(17) = pmd + dv(18) = pme + dv(19) = mfm + dv(20) = mfa + dv(21) = id2 + dv(22) = jmag + dv(23) = hmag + dv(24) = kmag + dv(25) = ejmag + dv(26) = ehmag + dv(27) = ekmag + dv(28) = iccj + dv(29) = icch + dv(30) = icck + dv(31) = phqj + dv(32) = phqh + dv(33) = phqk + dv(34) = apmag(1) + dv(35) = apmag(2) + dv(36) = apmag(3) + dv(37) = apmag(4) + dv(38) = apmag(5) + dv(39) = aperr(1) + dv(40) = aperr(2) + dv(41) = aperr(3) + dv(42) = aperr(4) + dv(43) = aperr(5) + dv(44) = no1 + dv(45) = mo1 + + RETURN ! normal return with valid data + + 91 eozf = .TRUE. ! end of file encountered + END ! subr. + +************************************************************************ + + SUBROUTINE xi_byte_flip (x0,zmax,bmax) +C +C input : x0 = array with index +C zmax= dimension of x0, max. number of zones +C bmax= dimension of x0, max. number of RA bins each zone +C output: x0 = same with byte flip applied + + INTEGER zmax,bmax + INTEGER x0 (zmax,bmax) + INTEGER zn,j + + DO zn= 1,zmax + DO j = 1,bmax + CALL flip4 (x0(zn,j)) + ENDDO + ENDDO + WRITE (*,'(a)') 'index array: byte flip done' + + END ! subr. + +************************************************************************ + + SUBROUTINE flip2 (i2) +C +C input: Integer*2 value i2 +C output: same with byte fliped + + IMPLICIT NONE + INTEGER*2 i2, i2i,i2o +C BYTE a(2), b(2) + INTEGER*1 a(2), b(2) + EQUIVALENCE (i2i,a) + EQUIVALENCE (i2o,b) + + i2i = i2 + b(1) = a(2) + b(2) = a(1) + i2 = i2o + + END ! subr. + +************************************************************************ + + SUBROUTINE flip4 (i4) +C +C input: Integer*4 value i4 +C output: same with byte fliped + + IMPLICIT NONE + INTEGER*4 i4, i4i,i4o +C BYTE a(4), b(4) + INTEGER*1 a(4), b(4) + EQUIVALENCE (i4i,a) + EQUIVALENCE (i4o,b) + + i4i = i4 + b(1) = a(4) + b(2) = a(3) + b(3) = a(2) + b(4) = a(1) + i4 = i4o + + END ! subr. + +************************************************************************ + + SUBROUTINE valid_range (dv,dmin,dmax) +C +C restrict given data value (dv) to given range (dmin,dmax) + + IMPLICIT NONE + REAL*8 dv, dmin,dmax + + IF (dv.LT.dmin) dv = dmin + IF (dv.GT.dmax) dv = dmax + + END ! subr. + +************************************************************************ + + SUBROUTINE get_zone_range (dc1,dc2,zmax, d1m,d2m,z1,z2,nz) +C +C input: dc1,dc2 = declination range (degree) +C zmax = largest zone number available +C output: d1m,d2m = declination range in mas +C z1, z2 = req. range of zone numbers (0.2 deg steps) +C nz = number of zones, or 0 if out of range + + IMPLICIT NONE + REAL*8 dc1,dc2, dcx + INTEGER zmax, d1m,d2m, z1,z2, nz + + IF (dc1.LT.-90.0d0.AND.dc2.LT.-90.0d0) THEN + nz = 0 + z1 = 1 + z2 = 0 + RETURN + ENDIF + + CALL valid_range (dc1,-90.0d0,90.0d0) + CALL valid_range (dc2,-90.0d0,90.0d0) + + IF (dc1.GT.dc2) THEN ! flip range + dcx = dc1 + dc1 = dc2 + dc2 = dcx + ENDIF + + d1m = IDNINT (dc1 * 3.6d6) ! declination degree to mas + d2m = IDNINT (dc2 * 3.6d6) + + z1 = (d1m + 324000000) / 720000 + 1 ! 720" = 0.2 deg + z2 = (d2m + 323999999) / 720000 + 1 + + IF (z2.GT.zmax) z2 = zmax + + IF (z1.GT.zmax) THEN ! out of available zone range + z1 = zmax + 1 + nz = 0 + ELSE + nz = z2 - z1 + 1 + ENDIF + + END ! subr. + +************************************************************************ + + SUBROUTINE get_ra_range (ra1,ra2, ralo,rahi,i1,i2,nr) +C +C input: ra1,ra2 = RA range (hour) +C output: ralo,rahi = range of RA in mas (1 or 2) +C i1, i2 = range in index for 1 arcmin bins +C nr = number of ranges = 1 or 2 +C 2 ranges possible, if ra1 > ra2, +C then assume cross over 24/0 hour in RA --> 2 ranges +C (e.g. 23...1 --> 23...24, and 0...1 hour for output) + + IMPLICIT NONE + REAL*8 ra1,ra2 + INTEGER ralo(2),rahi(2), i1(2),i2(2), nr + INTEGER r1m,r2m + REAL*8 rax + + CALL valid_range (ra1, 0.0d0,24.0d0) + CALL valid_range (ra2, 0.0d0,24.0d0) + + r1m = IDNINT (ra1 * 5.4d7) ! RA hour to mas + r2m = IDNINT (ra2 * 5.4d7) + + IF (r1m.LE.r2m) THEN ! normal case + nr = 1 + i1(1) = r1m / 900000 + 1 ! 900 arcsec = 0.25 deg + i2(1) = (r2m-1) / 900000 + 1 + i1(2) = 1 + i2(2) = 0 + ralo(1) = r1m + rahi(1) = r2m + ralo(2) = 0 + rahi(2) = 0 + + ELSE ! cross over 24/0 + nr = 2 + i1(1) = r1m / 900000 + 1 + i2(1) = 1440 ! = last bin up to 24 hour + i1(2) = 1 + i2(2) = (r2m-1) / 900000 + 1 + ralo(1) = r1m + rahi(1) = 1296000000 ! 24 hour in mas + ralo(2) = 0 + rahi(2) = r2m + ENDIF + + END ! subr. + +************************************************************************ + + SUBROUTINE as2hms (ra,dk,crekt,cdekl) +C +C convert R*8 RA, DC (arcsec) to hms, dms strings +C +C 940725 NZ update to CHARACTER*13 to 1/1000 arcsec + + IMPLICIT REAL*8 (A-H,L-Z) + + REAL*8 RA, DK + INTEGER*4 IRASTD, IRAMIN, IDKGRD,IDKMIN + CHARACTER*1 CVZ + CHARACTER*13 CREKT,CDEKL + + IF (RA.GT.1296000.D0) RA = RA - 1296000.D0 + IF (RA.LT. 0.D0) RA = RA + 1296000.D0 + IF (RA.GT.1296000.D0.OR.RA.LT.0.D0) THEN + WRITE (90,'(1X//1X,A,F13.3/)') 'RA > 24 hours or < 0 :',RA + RETURN + ENDIF + + RASTD = RA/(3600.D0*15.D0) + IRASTD= IDINT(RASTD) + RAREST= RASTD-DFLOAT(IRASTD) + IRAMIN= IDINT(RAREST*60.D0) + RASEC = RAREST*3600.D0-DFLOAT(IRAMIN)*60.D0 + + IF (DABS(RASEC-60.D0).LT.0.001D0) THEN + RASEC = 0.D0 + IRAMIN= IRAMIN+1 + IF (IRAMIN.EQ.60) THEN + IRAMIN= 0 + IRASTD= IRASTD+1 + IF (IRASTD.EQ.24) IRASTD= 0 + ENDIF + ENDIF + + IF (DABS(DK).GT.324000.D0) THEN + WRITE (90,'(1X//1X,A,F13.3/)') 'abs (DC) > 90 deg :',DK + RETURN + ENDIF + + DKGRD = DK/3600.D0 + CVZ= '+' + + IF (DK.LT.0.D0) THEN + CVZ= '-' + DKGRD= -DKGRD + ENDIF + + IDKGRD= IDINT(DKGRD) + DKREST= DKGRD-DFLOAT(IDKGRD) + IDKMIN= IDINT(DKREST*60.D0) + DKSEC = DKREST*3600.D0-DFLOAT(IDKMIN)*60.D0 + + IF (DABS(DKSEC-60.D0).LT.0.01D0) THEN + DKSEC = 0.D0 + IDKMIN= IDKMIN+1 + IF (IDKMIN.EQ.60) THEN + IDKMIN= 0 + IDKGRD= IDKGRD+1 + ENDIF + ENDIF + + WRITE (CREKT,'( I2.2,1X,I2.2,1X,F7.4)') + . IRASTD,IRAMIN,RASEC + IF (CREKT(7:7).EQ.' ') CREKT(7:7)= '0' + WRITE (CDEKL,'(A1,I2.2,1X,I2.2,1X,F6.3)') + . CVZ,IDKGRD,IDKMIN,DKSEC + IF (CDEKL(8:8).EQ.' ') CDEKL(8:8)= '0' + + RETURN + END + +************************************************************************ + + SUBROUTINE SORTI (A,IS,IW,IV,N,KX,KY) +C +C sort of 2-dimensional Integer array +C +C input: A = array +C IS = column index to sort by +C IW = number of columns to sort +C IV,N = lower, upper limit on line numbers to sort +C KX,KY= dimension of array A (lines,columns) +C output: A = same array after sort + + DIMENSION A(KX,KY),IL(60),IU(60),T(60),TT(60) + INTEGER A,T,TT + M=1 + I=IV + J=N + II=I + GOTO114 + 111 IJ=0.01+0.5*(I+J) + DO115IX=1,IW + 115 T(IX)=A(IJ,IX) + K=I + L=J + IF(A(I,IS).GT.T(IS))THEN + DO116IX=1,IW + A(IJ,IX)=A(I,IX) + A(I,IX)=T(IX) + 116 T(IX)=A(IJ,IX) + END IF + IF(A(J,IS).LT.T(IS))THEN + DO117IX=1,IW + A(IJ,IX)=A(J,IX) + A(J,IX)=T(IX) + 117 T(IX)=A(IJ,IX) + IF(A(I,IS).GT.T(IS))THEN + DO118IX=1,IW + A(IJ,IX)=A(I,IX) + A(I,IX)=T(IX) + 118 T(IX)=A(IJ,IX) + END IF + END IF + 112 L=L-1 + IF(A(L,IS).GT.T(IS))GOTO112 + DO119IX=1,IW + 119 TT(IX)=A(L,IX) + 113 K=K+1 + IF(A(K,IS).LT.T(IS))GOTO113 + IF(K.LE.L)THEN + DO120IX=1,IW + A(L,IX)=A(K,IX) + 120 A(K,IX)=TT(IX) + GOTO112 + END IF + IF((L-I).GT.(J-K))THEN + IL(M)=I + IU(M)=L + I=K + ELSE + IL(M)=K + IU(M)=J + J=L + END IF + M=M+1 + 114 IF(J-I.GT.10)GOTO 111 + IF(I.EQ.II)THEN + IF(I.LT.J)GOTO111 + END IF + NI=I+1 + DO121IZ=NI,J + I=IZ + DO122IX=1,IW + 122 T(IX)=A(IZ,IX) + K=I-1 + IF(A(K,IS).GT.T(IS))THEN + 123 DO124IX=1,IW + 124 A(K+1,IX)=A(K,IX) + K=K-1 + IF(A(K,IS).GT.T(IS))GOTO123 + DO125IX=1,IW + 125 A(K+1,IX)=T(IX) + END IF + 121 CONTINUE + M=M-1 + IF(M.GE.1)THEN + I=IL(M) + J=IU(M) + GOTO114 + END IF + RETURN + END ! subr. + From f0375838f88aa9ee5b1068e5d68dfdf3c2d6631a Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Thu, 25 Jul 2019 14:05:42 -0700 Subject: [PATCH 03/30] mild update to fortran binaries --- py/desitarget/fortran/ADM-v1dump.f | 2 +- py/desitarget/fortran/v1dump | Bin 31064 -> 31064 bytes 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/py/desitarget/fortran/ADM-v1dump.f b/py/desitarget/fortran/ADM-v1dump.f index 751ec28a7..8a0d7712a 100644 --- a/py/desitarget/fortran/ADM-v1dump.f +++ b/py/desitarget/fortran/ADM-v1dump.f @@ -38,7 +38,7 @@ PROGRAM v1dump patho = TRIM(TRIM(uratdir) // "/csv/") C ADM the file RANGE to convert. zn1 = 326 - zn2 = 327 + zn2 = 900 C ADM byte-flip or not (big or little endian). bfc = 'N' C ADM the delimiter for the output files. diff --git a/py/desitarget/fortran/v1dump b/py/desitarget/fortran/v1dump index 025d0da9d2874b0a0cc8b9d6015f7f93d836383e..4ea0d2887681ca066e88ebe2ddc42942028a2c94 100755 GIT binary patch delta 59 zcmccdiSfoK#tqv9nOc}P?+~2ipkkP8kYr+0|lD0W7=#?4(ZMqX- From eeb1652555f076f9fc7497e717f1a04589e3484b Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Thu, 25 Jul 2019 14:07:29 -0700 Subject: [PATCH 04/30] more minor updates to fortran --- py/desitarget/fortran/ADM-v1dump.f | 7 +++---- py/desitarget/fortran/v1dump | Bin 31064 -> 31064 bytes 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/py/desitarget/fortran/ADM-v1dump.f b/py/desitarget/fortran/ADM-v1dump.f index 8a0d7712a..4d35c7dfb 100644 --- a/py/desitarget/fortran/ADM-v1dump.f +++ b/py/desitarget/fortran/ADM-v1dump.f @@ -1,9 +1,9 @@ - PROGRAM v1dump + PROGRAM v1dump C C gfortran -o v1dump ADM-v1dump.f v1sub.f C C - loop over range of zones -C - read individual URAT zone file (binary data) +C - read individual URAT zone file (binary data) C - output all data to ASCII file, no header lines C - user select column separator character, no empty columns C - for explanations of the 45 columns see "readme" file @@ -15,7 +15,7 @@ PROGRAM v1dump C 100725 ADM Adapted for NERSC and hardcoded to skip interactive steps IMPLICIT NONE INTEGER dimc, dims - PARAMETER (dimc = 45 ! number of URAT data columns + PARAMETER (dimc = 45 ! number of URAT data columns . ,dims = 45) ! number of separators INTEGER dv(dimc), mi(dimc), ma(dimc), is(dims) @@ -94,4 +94,3 @@ PROGRAM v1dump ENDDO ! loop all zones END ! main - diff --git a/py/desitarget/fortran/v1dump b/py/desitarget/fortran/v1dump index 4ea0d2887681ca066e88ebe2ddc42942028a2c94..389e8e60a564de90cf64716480ee21689e31a92f 100755 GIT binary patch delta 51 zcmccdiSfoK#tmLUDi$fmhH1&EmbyvFY38~nh6Wb8Nd~4Cy2)wDDdwiBDM^V2shjhH HUa0{9=2Q{z delta 51 zcmccdiSfoK#tmLUDu&4hNhZc=rn*T+hNikEhDjE>7AEG#y2(i Date: Tue, 30 Jul 2019 15:04:19 -0700 Subject: [PATCH 05/30] switch to reading healpix files when calculating pixweight stellar densities, to limit I/O --- py/desitarget/randoms.py | 35 +++++++++++++++++++++++------------ 1 file changed, 23 insertions(+), 12 deletions(-) diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py index b2d45d5dd..67a213057 100644 --- a/py/desitarget/randoms.py +++ b/py/desitarget/randoms.py @@ -741,28 +741,39 @@ def stellar_density(nside=256): """ # ADM check that the GAIA_DIR is set and retrieve it. gaiadir = _get_gaia_dir() - fitsdir = os.path.join(gaiadir, 'fits') + hpdir = os.path.join(gaiadir, 'healpix') - # ADM the number of pixels and the pixel area at the passed nside. + # ADM the number of pixels and the pixel area at nside. npix = hp.nside2npix(nside) pixarea = hp.nside2pixarea(nside, degrees=True) - # ADM an output array to populate containing all possible HEALPixels at the passed nside. + # ADM an output array of all possible HEALPixels at nside. pixout = np.zeros(npix, dtype='int32') # ADM find all of the Gaia files. - filenames = glob(fitsdir+'/*fits') + filenames = glob(os.path.join(hpdir,'*fits')) - # ADM read in each file, restricting to the criteria for point sources - # ADM and storing in a HEALPixel map at resolution nside. - for filename in filenames: - # ADM save memory and speed up by only reading in a subset of columns. - gobjs = fitsio.read(filename, - columns=['RA', 'DEC', 'PHOT_G_MEAN_MAG', 'ASTROMETRIC_EXCESS_NOISE']) + # ADM read in each file, restricting to the criteria for point + # ADM sources and storing in a HEALPixel map at resolution nside. + nfiles = len(filenames) + t0 = time() + for nfile, filename in enumerate(filenames): + if nfile % 200 == 0 and nfile > 0: + elapsed = time() - t0 + rate = elapsed / nfile + log.info('{}/{} files; {:.1f} secs/file; {:.1f} total mins elapsed' + .format(nfile, nfiles, rate, elapsed/60.)) + + # ADM save memory, speed up by only reading a subset of columns. + gobjs = fitsio.read( + filename, + columns=['RA', 'DEC', 'PHOT_G_MEAN_MAG', 'ASTROMETRIC_EXCESS_NOISE'] + ) - # ADM restrict to subset of sources using point source definition. + # ADM restrict to subset of point sources. ra, dec = gobjs["RA"], gobjs["DEC"] - gmag, excess = gobjs["PHOT_G_MEAN_MAG"], gobjs["ASTROMETRIC_EXCESS_NOISE"] + gmag = gobjs["PHOT_G_MEAN_MAG"] + excess = gobjs["ASTROMETRIC_EXCESS_NOISE"] point = (excess == 0.) | (np.log10(excess) < 0.3*gmag-5.3) grange = (gmag >= 12) & (gmag < 17) w = np.where(point & grange) From e55127173ebb8048123dbc01b0c7eb28e3d9c5e2 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Tue, 30 Jul 2019 15:26:23 -0700 Subject: [PATCH 06/30] default to creating an 'all-sky' GFAs file --- bin/select_gfas | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/bin/select_gfas b/bin/select_gfas index 6117fa3ee..9067bc81f 100755 --- a/bin/select_gfas +++ b/bin/select_gfas @@ -32,7 +32,7 @@ ap.add_argument('-s2', "--surveydir2", help='Additional Legacy Surveys Data Release directory (useful for combining, e.g., DR8 into one file of GFAs)', default=None) ap.add_argument('-t', "--tiles", - help="File specifying the tiles to which to restrict the GFAs (defaults to all tiles in the DESI footprint)", + help="File specifying the tiles to which to restrict the GFAs", default=None) ap.add_argument('-dec', "--mindec", type=float, help="Minimum declination to include in output file for NON-LEGACY-SURVEYS sources (degrees; defaults to [-30])", @@ -40,15 +40,23 @@ ap.add_argument('-dec', "--mindec", type=float, ap.add_argument('-b', "--mingalb", type=float, help="Closest latitude to Galactic plane to output for NON-LEGACY-SURVEYS sources (e.g. send the default [10] to limit to areas beyond -10o <= b < 10o)", default=10) -ap.add_argument("--cmx", action='store_true', - help="If sent, then create a commissioning file. Commissioning files are not restricted by any tiling pattern, even if --tiles is sent") +ap.add_argument("--desifootprint", action='store_true', + help="If sent, then limit to the current DESIMODEL Main Survey footprint (the default is to produce an 'all-sky' file)") ns = ap.parse_args() -survey = 'main' -if ns.cmx: - log.info('Making a cmx-style file. No restrictions on tiling') +if ns.tiles is not None: + log.info('Restricting to passed tiles') + survey = 'tiles' + cmx = False +elif ns.desifootprint: + log.info('Restricting to the DESI Main Survey footprint') + survey = 'main' + cmx = False +else: + log.info('Producing GFAs across the entire sky (with specified limits)') survey = 'cmx' + cmx = True infiles = io.list_sweepfiles(ns.surveydir) if ns.surveydir2 is not None: @@ -66,7 +74,7 @@ if len(infiles) == 0: log.info('running on {} processors...t = {:.1f} mins'.format(ns.numproc, (time()-t0)/60.)) gfas = select_gfas(infiles, maglim=ns.maglim, numproc=ns.numproc, tilesfile=ns.tiles, - cmx=ns.cmx, mindec=ns.mindec, mingalb=ns.mingalb) + cmx=cmx, mindec=ns.mindec, mingalb=ns.mingalb) io.write_gfas(ns.dest, gfas, indir=ns.surveydir, indir2=ns.surveydir2, nside=nside, survey=survey) From ff24582f8c2a1f52e25efc769483c83fdad53f70 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Tue, 30 Jul 2019 18:06:41 -0700 Subject: [PATCH 07/30] Add a README file explaining the fortran directory --- py/desitarget/fortran/README | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 py/desitarget/fortran/README diff --git a/py/desitarget/fortran/README b/py/desitarget/fortran/README new file mode 100644 index 000000000..9d0bafd2b --- /dev/null +++ b/py/desitarget/fortran/README @@ -0,0 +1,17 @@ +The files in this directory are: + +ADM-v1dump.f: FORTRAN code to parse UCAC/URAT files, adapted from + http://cdsarc.u-strasbg.fr/ftp/I/329/URAT1/access/v1dump.f + by Adam D. Myers (University of Wyoming). + +v1sub.f: Routines to help parse UCAC/URAT files, taken directly from + http://cdsarc.u-strasbg.fr/ftp/I/329/URAT1/access/v1sub.f + +v1dump: An executable compiled at NERSC via: + + gfortran -o v1dump ADM-v1dump.f v1sub.f + + +More about the UCAC/URAT survey is available at: + +http://cdsarc.u-strasbg.fr/viz-bin/Cat?I/329 From 8e21158998db7009f2f5fb8b6b728f213fb3d28f Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Wed, 31 Jul 2019 09:09:36 -0700 Subject: [PATCH 08/30] a little more logging when generating pixweight files --- py/desitarget/randoms.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py index 67a213057..0ae0161fa 100644 --- a/py/desitarget/randoms.py +++ b/py/desitarget/randoms.py @@ -760,8 +760,8 @@ def stellar_density(nside=256): for nfile, filename in enumerate(filenames): if nfile % 200 == 0 and nfile > 0: elapsed = time() - t0 - rate = elapsed / nfile - log.info('{}/{} files; {:.1f} secs/file; {:.1f} total mins elapsed' + rate = nfile / elapsed + log.info('{}/{} files; {:.1f} files/sec; {:.1f} total mins elapsed' .format(nfile, nfiles, rate, elapsed/60.)) # ADM save memory, speed up by only reading a subset of columns. @@ -978,18 +978,26 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None): pixels, pixcnts = np.unique(pixnums, return_counts=True) pixcnts = np.insert(pixcnts, 0, 0) pixcnts = np.cumsum(pixcnts) - + log.info('Done sorting...t = {:.1f}s'.format(time()-start)) # ADM work through the ordered pixels to populate the median for # ADM each quantity of interest. cols = ['EBV', 'PSFDEPTH_W1', 'PSFDEPTH_W2', 'PSFDEPTH_G', 'GALDEPTH_G', 'PSFSIZE_G', 'PSFDEPTH_R', 'GALDEPTH_R', 'PSFSIZE_R', 'PSFDEPTH_Z', 'GALDEPTH_Z', 'PSFSIZE_Z'] - for i in range(len(pixcnts)-1): + t0 = time() + npix = len(pixcnts) + stepper = npix//50 + for i in range(npix-1): inds = pixorder[pixcnts[i]:pixcnts[i+1]] pix = pixnums[inds][0] for col in cols: hpxinfo[col][pix] = np.median(randoms[col][inds]) + if i % stepper == 0 and i > 0: + elapsed = time() - t0 + rate = i / elapsed + log.info('{}/{} pixels; {:.1f} pix/sec; {:.1f} total mins elapsed' + .format(i, npix, rate, elapsed/60.)) log.info('Done...t = {:.1f}s'.format(time()-start)) From bf233c554db1c3160d011fe693e115e182694d22 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Wed, 31 Jul 2019 16:11:25 -0700 Subject: [PATCH 09/30] update RA/Dec matcher to return the match separation, if requested --- py/desitarget/gaiamatch.py | 2 +- py/desitarget/geomask.py | 11 ++++++++++- py/desitarget/randoms.py | 2 +- 3 files changed, 12 insertions(+), 3 deletions(-) diff --git a/py/desitarget/gaiamatch.py b/py/desitarget/gaiamatch.py index 765cb1ef3..00f72c350 100644 --- a/py/desitarget/gaiamatch.py +++ b/py/desitarget/gaiamatch.py @@ -770,7 +770,7 @@ def match_gaia_to_primary(objs, matchrad=1., retaingaia=False, ------- :class:`~numpy.ndarray` The matching Gaia information for each object, where the returned format and - columns correspond to `desitarget.secondary.gaiadatamodel` + columns correspond to `desitarget.gaiamatch.gaiadatamodel` Notes ----- diff --git a/py/desitarget/geomask.py b/py/desitarget/geomask.py index 620cda5fe..74312ee78 100644 --- a/py/desitarget/geomask.py +++ b/py/desitarget/geomask.py @@ -1256,7 +1256,7 @@ def is_in_gal_box(objs, lbbox, radec=False): return ii -def radec_match_to(matchto, objs, sep=1., radec=False): +def radec_match_to(matchto, objs, sep=1., radec=False, return_sep=False): """Match objects to a catalog list on RA/Dec. Parameters @@ -1270,6 +1270,9 @@ def radec_match_to(matchto, objs, sep=1., radec=False): radec : :class:`bool`, optional, defaults to ``False`` If ``True`` then `objs` and `matchto` are [RA, Dec] lists instead of rec arrays. + return_sep : :class:`bool`, optional, defaults to ``False`` + If ``True`` then return the separation between each object, not + just the indexes of the match. Returns ------- @@ -1277,6 +1280,9 @@ def radec_match_to(matchto, objs, sep=1., radec=False): The indexes in `match2` for which `objs` matches `match2` at < `sep`. :class:`~numpy.ndarray` (of integers) The indexes in `objs` for which `objs` matches `match2` at < `sep`. + :class:`~numpy.ndarray` (of floats) + The distances in ARCSECONDS of the matches. + Only returned if `return_sep` is ``True``. Notes ----- @@ -1313,4 +1319,7 @@ def radec_match_to(matchto, objs, sep=1., radec=False): ii = d2d < sep*u.arcsec + if return_sep: + return idmatchto[ii], idobjs[ii], d2d[ii].arcsec + return idmatchto[ii], idobjs[ii] diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py index 0ae0161fa..faa9b7e05 100644 --- a/py/desitarget/randoms.py +++ b/py/desitarget/randoms.py @@ -758,7 +758,7 @@ def stellar_density(nside=256): nfiles = len(filenames) t0 = time() for nfile, filename in enumerate(filenames): - if nfile % 200 == 0 and nfile > 0: + if nfile % 1000 == 0 and nfile > 0: elapsed = time() - t0 rate = nfile / elapsed log.info('{}/{} files; {:.1f} files/sec; {:.1f} total mins elapsed' From 7bc8dd6f4c184ed6724ec62a2667bd4cbe9b9e1c Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Wed, 31 Jul 2019 16:12:01 -0700 Subject: [PATCH 10/30] new module to download and format URAT files, and match to URAT locations --- py/desitarget/uratmatch.py | 603 +++++++++++++++++++++++++++++++++++++ 1 file changed, 603 insertions(+) create mode 100644 py/desitarget/uratmatch.py diff --git a/py/desitarget/uratmatch.py b/py/desitarget/uratmatch.py new file mode 100644 index 000000000..90f273776 --- /dev/null +++ b/py/desitarget/uratmatch.py @@ -0,0 +1,603 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +# -*- coding: utf-8 -*- +""" +==================== +desitarget.uratmatch +==================== + +Useful `URAT`_ matching and manipulation routines. + +.. _`URAT`: http://cdsarc.u-strasbg.fr/viz-bin/cat/I/329 +""" +import os +import numpy as np +import fitsio +import requests +import pickle + +from pkg_resources import resource_filename +from time import time +from astropy.io import ascii +from glob import glob +import healpy as hp + +from desitarget.internal import sharedmem +from desimodel.footprint import radec2pix +from desitarget.geomask import add_hp_neighbors, radec_match_to + +# ADM set up the DESI default logger +from desiutil.log import get_logger +log = get_logger() + +# ADM start the clock +start = time() + +# ADM columns contained in our version of the URAT fits files. +uratdatamodel = np.array([], dtype=[ + ('URAT_ID', '>i8'), ('RA', '>f8'), ('DEC', '>f8'), + ('APASS_G_MAG', '>f4'), ('APASS_G_MAG_ERROR', '>f4'), + ('APASS_R_MAG', '>f4'), ('APASS_R_MAG_ERROR', '>f4'), + ('APASS_I_MAG', '>f4'), ('APASS_I_MAG_ERROR', '>f4'), + ('PMRA', '>f4'), ('PMDEC', '>f4'), ('PM_ERROR', '>f4') +]) + +def _get_urat_dir(): + """Convenience function to grab the URAT environment variable. + + Returns + ------- + :class:`str` + The directory stored in the $URAT_DIR environment variable. + """ + # ADM check that the $URAT_DIR environment variable is set. + uratdir = os.environ.get('URAT_DIR') + if uratdir is None: + msg = "Set $URAT_DIR environment variable!" + log.critical(msg) + raise ValueError(msg) + + return uratdir + + +def _get_urat_nside(): + """Grab the HEALPixel nside to be used throughout this module. + + Returns + ------- + :class:`int` + The HEALPixel nside number for URAT file creation and retrieval. + """ + nside = 32 + + return nside + + +def scrape_urat(url="http://cdsarc.u-strasbg.fr/ftp/I/329/URAT1/v12/", + nfiletest=None): + """Retrieve the binary versions of the URAT files. + + Parameters + ---------- + url : :class:`str` + The web directory that hosts the archived binary URAT files. + nfiletest : :class:`int`, optional, defaults to ``None`` + If an integer is sent, only retrieve this number of files, for testing. + + Returns + ------- + Nothing + But the archived URAT files are written to $URAT_DIR/binary. + + Notes + ----- + - The environment variable $URAT_DIR must be set. + - Runs in about 50 minutes for 575 URAT files. + """ + # ADM check that the URAT_DIR is set and retrieve it. + uratdir = _get_urat_dir() + + # ADM construct the directory to which to write files. + bindir = os.path.join(uratdir, 'binary') + # ADM the directory better be empty for the wget! + if os.path.exists(bindir): + if len(os.listdir(bindir)) > 0: + msg = "{} should be empty to wget URAT binary files!".format(bindir) + log.critical(msg) + raise ValueError(msg) + # ADM make the directory, if needed. + else: + log.info('Making URAT directory for storing binary files') + os.makedirs(bindir) + + index = requests.get(url) + + # ADM retrieve any file name that starts with z. + # ADM the [1::2] pulls back just the odd lines from the split list. + garbled = index.text.split("z")[1::2] + filelist = ["z{}".format(g[:3]) for g in garbled] + + # ADM if nfiletest was passed, just work with that number of files. + test = nfiletest is not None + if test: + filelist = filelist[:nfiletest] + nfiles = len(filelist) + + # ADM loop through the filelist. + start = time() + for nfile, fileinfo in enumerate(filelist): + # ADM make the wget command to retrieve the file and issue it. + cmd = 'wget -q {} -P {}'.format(os.path.join(url, fileinfo), bindir) + print(cmd) + os.system(cmd) + if nfile % 25 == 0 or test: + elapsed = time() - start + rate = nfile / elapsed + log.info( + '{}/{} files; {:.1f} files/sec; {:.1f} total mins elapsed' + .format(nfile+1, nfiles, rate, elapsed/60.) + ) + + log.info('Done...t={:.1f}s'.format(time()-start)) + + return + + +def urat_binary_to_csv(): + """Convert files in $URAT_DIR/binary to files in $URAT_DIR/csv. + + Returns + ------- + Nothing + But the archived URAT binary files in $URAT_DIR/binary are + converted to CSV files in the $URAT_DIR/csv. + + Notes + ----- + - The environment variable $URAT_DIR must be set. + - Relies on the executable fortran/v1dump, which is only + tested at NERSC. + - Runs in about 40 minutes for 575 files. + """ + # ADM check that the URAT_DIR is set. + uratdir = _get_urat_dir() + + log.info('Begin converting URAT files to CSV...t={:.1f}s' + .format(time()-start)) + + cmd = resource_filename('desitarget', 'fortran/v1dump') + os.system(cmd) + + log.info('Done...t={:.1f}s'.format(time()-start)) + + return + + +def urat_csv_to_fits(numproc=5): + """Convert files in $URAT_DIR/csv to files in $URAT_DIR/fits. + + Parameters + ---------- + numproc : :class:`int`, optional, defaults to 5 + The number of parallel processes to use. + + Returns + ------- + Nothing + But the archived URAT CSV files in $URAT_DIR/csv are converted + to FITS files in the directory $URAT_DIR/fits. Also, a look-up + table is written to $URAT_DIR/fits/hpx-to-files.pickle for which + each index is an nside=_get_urat_nside(), nested scheme HEALPixel + and each entry is a list of the FITS files that touch that HEAPixel. + + Notes + ----- + - The environment variable $URAT_DIR must be set. + - if numproc==1, use the serial code instead of the parallel code. + - Runs in about 10 minutes with numproc=25 for 575 files. + """ + # ADM the resolution at which the URAT HEALPix files should be stored. + nside = _get_urat_nside() + + # ADM check that the URAT_DIR is set. + uratdir = _get_urat_dir() + log.info("running on {} processors".format(numproc)) + + # ADM construct the directories for reading/writing files. + csvdir = os.path.join(uratdir, 'csv') + fitsdir = os.path.join(uratdir, 'fits') + + # ADM make sure the output directory is empty. + if os.path.exists(fitsdir): + if len(os.listdir(fitsdir)) > 0: + msg = "{} should be empty to make URAT FITS files!".format(fitsdir) + log.critical(msg) + raise ValueError(msg) + # ADM make the output directory, if needed. + else: + log.info('Making URAT directory for storing FITS files') + os.makedirs(fitsdir) + + # ADM construct the list of input files. + infiles = glob("{}/*csv*".format(csvdir)) + nfiles = len(infiles) + + # ADM the critical function to run on every file. + def _write_urat_fits(infile): + """read an input name for a csv file and write it to FITS""" + outbase = os.path.basename(infile) + outfilename = "{}.fits".format(outbase.split(".")[0]) + outfile = os.path.join(fitsdir, outfilename) + # ADM astropy understands without specifying format='csv'. + fitstable = ascii.read(infile) + + # ADM map the ascii-read csv to typical DESI quantities. + nobjs = len(fitstable) + done = np.zeros(nobjs, dtype=uratdatamodel.dtype) + # ADM have to do this one-by-one, given the format. + done["RA"] = fitstable['col1']/1000./3600. + done["DEC"] = fitstable['col2']/1000./3600. - 90. + done["PMRA"] = fitstable['col16']/10. + done["PMDEC"] = fitstable['col17']/10. + done["PM_ERROR"] = fitstable['col18']/10. + done["APASS_G_MAG"] = fitstable['col36']/1000. + done["APASS_R_MAG"] = fitstable['col37']/1000. + done["APASS_I_MAG"] = fitstable['col38']/1000. + done["APASS_G_MAG_ERROR"] = fitstable['col41']/1000. + done["APASS_R_MAG_ERROR"] = fitstable['col42']/1000. + done["APASS_I_MAG_ERROR"] = fitstable['col43']/1000. + done["URAT_ID"] = fitstable['col46'] + + fitsio.write(outfile, done, extname='URATFITS') + + # ADM return the HEALPixels that this file touches. + pix = set(radec2pix(nside, done["RA"], done["DEC"])) + return [pix, os.path.basename(outfile)] + + # ADM this is just to count processed files in _update_status. + nfile = np.zeros((), dtype='i8') + t0 = time() + + def _update_status(result): + """wrapper function for the critical reduction operation, + that occurs on the main parallel process""" + if nfile % 25 == 0 and nfile > 0: + rate = nfile / (time() - t0) + elapsed = time() - t0 + log.info( + '{}/{} files; {:.1f} files/sec; {:.1f} total mins elapsed' + .format(nfile, nfiles, rate, elapsed/60.) + ) + nfile[...] += 1 # this is an in-place modification + return result + + # - Parallel process input files... + if numproc > 1: + pool = sharedmem.MapReduce(np=numproc) + with pool: + pixinfile = pool.map(_write_urat_fits, infiles, reduce=_update_status) + # ADM ...or run in serial. + else: + pixinfile = list() + for file in infiles: + pixinfile.append(_update_status(_write_urat_fits(file))) + + # ADM create a list for which each index is a HEALPixel and each + # ADM entry is a list of files that touch that HEALPixel. + npix = hp.nside2npix(nside) + pixlist = [[] for i in range(npix)] + for pixels, file in pixinfile: + for pix in pixels: + pixlist[pix].append(file) + + # ADM write out the HEALPixel->files look-up table. + outfilename = os.path.join(fitsdir, "hpx-to-files.pickle") + outfile = open(outfilename, "wb") + pickle.dump(pixlist, outfile) + outfile.close() + + log.info('Done...t={:.1f}s'.format(time()-t0)) + + return + + +def urat_fits_to_healpix(numproc=5): + """Convert files in $URAT_DIR/fits to files in $URAT_DIR/healpix. + + Parameters + ---------- + numproc : :class:`int`, optional, defaults to 5 + The number of parallel processes to use. + + Returns + ------- + Nothing + But the archived URAT FITS files in $URAT_DIR/fits are + rearranged by HEALPixel in the directory $URAT_DIR/healpix. + The HEALPixel sense is nested with nside=_get_urat_nside(), and + each file in $URAT_DIR/healpix is called healpix-xxxxx.fits, + where xxxxx corresponds to the HEALPixel number. + + Notes + ----- + - The environment variable $URAT_DIR must be set. + - if numproc==1, use the serial code instead of the parallel code. + - Runs in about 10 minutes with numproc=25. + """ + # ADM the resolution at which the URAT HEALPix files should be stored. + nside = _get_urat_nside() + + # ADM check that the URAT_DIR is set. + uratdir = _get_urat_dir() + + # ADM construct the directories for reading/writing files. + fitsdir = os.path.join(uratdir, 'fits') + hpxdir = os.path.join(uratdir, 'healpix') + + # ADM make sure the output directory is empty. + if os.path.exists(hpxdir): + if len(os.listdir(hpxdir)) > 0: + msg = "{} should be empty to make URAT HEALPix files!".format(hpxdir) + log.critical(msg) + raise ValueError(msg) + # ADM make the output directory, if needed. + else: + log.info('Making URAT directory for storing HEALPix files') + os.makedirs(hpxdir) + + # ADM read the pixel -> file look-up table. + infilename = os.path.join(fitsdir, "hpx-to-files.pickle") + infile = open(infilename, "rb") + pixlist = pickle.load(infile) + npixels = len(pixlist) + # ADM include the pixel number explicitly in the look-up table. + pixlist = list(zip(np.arange(npixels), pixlist)) + + # ADM the critical function to run on every file. + def _write_hpx_fits(pixlist): + """from files that touch a pixel, write out objects in each pixel""" + pixnum, files = pixlist + # ADM only proceed if some files touch a pixel. + if len(files) > 0: + # ADM track if it's our first time through the files loop. + first = True + # ADM Read in files that touch a pixel. + for file in files: + filename = os.path.join(fitsdir, file) + objs = fitsio.read(filename) + # ADM only retain objects in the correct pixel. + pix = radec2pix(nside, objs["RA"], objs["DEC"]) + if first: + done = objs[pix == pixnum] + first = False + else: + done = np.hstack([done, objs[pix == pixnum]]) + # ADM construct the name of the output file. + outfilename = 'healpix-{:05d}.fits'.format(pixnum) + outfile = os.path.join(hpxdir, outfilename) + # ADM write out the file. + hdr = fitsio.FITSHDR() + hdr['HPXNSIDE'] = nside + hdr['HPXNEST'] = True + fitsio.write(outfile, done, extname='URATHPX', header=hdr) + + return + + # ADM this is just to count processed files in _update_status. + npix = np.zeros((), dtype='i8') + t0 = time() + + def _update_status(result): + """wrapper function for the critical reduction operation, + that occurs on the main parallel process""" + if npix % 500 == 0 and npix > 0: + rate = npix / (time() - t0) + elapsed = time() - t0 + log.info( + '{}/{} files; {:.1f} files/sec; {:.1f} total mins elapsed' + .format(npix, npixels, rate, elapsed/60.) + ) + npix[...] += 1 # this is an in-place modification + return result + + # - Parallel process input files... + if numproc > 1: + pool = sharedmem.MapReduce(np=numproc) + with pool: + _ = pool.map(_write_hpx_fits, pixlist, reduce=_update_status) + # ADM ...or run in serial. + else: + for pix in pixlist: + _update_status(_write_hpx_fits(pix)) + + log.info('Done...t={:.1f}s'.format(time()-t0)) + + return + + +def make_urat_files(numproc=5, download=False): + """Make the HEALPix-split URAT files in one fell swoop. + + Parameters + ---------- + numproc : :class:`int`, optional, defaults to 5 + The number of parallel processes to use. + download : :class:`bool`, optional, defaults to ``False`` + If ``True`` then wget the URAT binary files from Vizier. + + Returns + ------- + Nothing + But produces: + - URAT DR1 binary files in $URAT_DIR/binary (if download=True). + - URAT CSV files with all URAT columns in $URAT_DIR/csv. + - FITS files with columns from `uratdatamodel` in $URAT_DIR/fits. + - FITS files reorganized by HEALPixel in $URAT_DIR/healpix. + + The HEALPixel sense is nested with nside=_get_urat_nside(), and + each file in $URAT_DIR/healpix is called healpix-xxxxx.fits, + where xxxxx corresponds to the HEALPixel number. + + Notes + ----- + - The environment variable $URAT_DIR must be set. + - if numproc==1, use the serial, instead of the parallel, code. + - Runs in about 2 hours with numproc=25 if download is ``True``. + - Runs in about 1 hour with numproc=25 if download is ``False``. + """ + t0 = time() + log.info('Begin making URAT files...t={:.1f}s'.format(time()-t0)) + + # ADM check that the URAT_DIR is set. + uratdir = _get_urat_dir() + + # ADM a quick check that the fits and healpix directories are empty + # ADM before embarking on the slower parts of the code. + csvdir = os.path.join(uratdir, 'csv') + fitsdir = os.path.join(uratdir, 'fits') + hpxdir = os.path.join(uratdir, 'healpix') + for direc in [csvdir, fitsdir, hpxdir]: + if os.path.exists(direc): + if len(os.listdir(direc)) > 0: + msg = "{} should be empty to make URAT files!".format(direc) + log.critical(msg) + raise ValueError(msg) + + if download: + scrape_urat() + log.info('Retrieved URAT files from Vizier...t={:.1f}s' + .format(time()-t0)) + + urat_binary_to_csv(numproc=numproc) + log.info('Converted binary files to CSV...t={:.1f}s'.format(time()-t0)) + + urat_csv_to_fits(numproc=numproc) + log.info('Converted CSV files to FITS...t={:.1f}s'.format(time()-t0)) + + urat_fits_to_healpix(numproc=numproc) + log.info('Rearranged FITS files by HEALPixel...t={:.1f}s'.format(time()-t0)) + + return + + +def find_urat_files(objs, neighbors=True, radec=False): + """Find full paths to URAT healpix files for objects by RA/Dec. + + Parameters + ---------- + objs : :class:`~numpy.ndarray` + Array of objects. Must contain the columns "RA" and "DEC". + neighbors : :class:`bool`, optional, defaults to ``True`` + Also return all pixels that touch the files of interest + to prevent edge effects (e.g. if a URAT source is 1 arcsec + away from a primary source and so in an adjacent pixel). + radec : :class:`bool`, optional, defaults to ``False`` + If ``True`` then the passed `objs` is an [RA, Dec] list + instead of a rec array that contains "RA" and "DEC". + + Returns + ------- + :class:`list` + A list of all URAT files to read to account for objects at + the passed locations. + + Notes + ----- + - The environment variable $URAT_DIR must be set. + """ + # ADM the resolution at which the URAT HEALPix files are stored. + nside = _get_urat_nside() + + # ADM check that the URAT_DIR is set and retrieve it. + uratdir = _get_urat_dir() + hpxdir = os.path.join(uratdir, 'healpix') + + # ADM which flavor of RA/Dec was passed. + if radec: + ra, dec = objs + else: + ra, dec = objs["RA"], objs["DEC"] + + # ADM convert RA/Dec to co-latitude and longitude in radians. + theta, phi = np.radians(90-dec), np.radians(ra) + + # ADM retrieve the pixels in which the locations lie. + pixnum = hp.ang2pix(nside, theta, phi, nest=True) + # ADM if neighbors was sent, then retrieve all pixels that touch each + # ADM pixel covered by the provided locations, to prevent edge effects... + if neighbors: + pixnum = add_hp_neighbors(nside, pixnum) + + # ADM reformat file names in the URAT healpix format. + uratfiles = [os.path.join(hpxdir, 'healpix-{:05d}.fits'.format(pn)) + for pn in pixnum] + + # ADM restrict to only files/HEALPixels actually covered by URAT. + uratfiles = [fn for fn in uratfiles if os.path.exists(fn)] + + return uratfiles + + +def match_to_urat(objs, matchrad=1.): + """Match objects to URAT healpix files and return URAT information. + + Parameters + ---------- + objs : :class:`~numpy.ndarray` + Must contain at least "RA" and "DEC". + matchrad : :class:`float`, optional, defaults to 1 arcsec + The radius at which to match in arcseconds. + + Returns + ------- + :class:`~numpy.ndarray` + The matching URAT information for each object. The returned + format is as for desitarget.uratmatch.uratdatamodel with + and extra column "DISTANCE" which is the matching distance + in ARCSECONDS. + + Notes + ----- + - For objects that do NOT have a match in URAT, the "URAT_ID" + and "DISTANCE" columns are -1, and other columns are zero. + - Retrieves the CLOSEST match to URAT for each passed object. + - Because this reads in HEALPixel split files, it's (far) faster + for objects that are clumped rather than widely distributed. + """ + # ADM set up an array of URAT information for the output. + nobjs = len(objs) + done = np.zeros(nobjs, dtype=uratdatamodel.dtype) + + # ADM objects without matches should have URAT_ID, DISTANCE of -1. + done["URAT_ID"] = -1 + distance = np.zeros(nobjs) - 1 + + # ADM determine which URAT files need to be scraped. + uratfiles = find_urat_files(objs) + nfiles = len(uratfiles) + + # ADM catch the case of no matches to URAT. + if nfiles == 0: + return done + + # ADM loop through the URAT files and find matches. + for ifn, fn in enumerate(uratfiles): + if ifn % 100 == 0 and ifn > 0: + log.info('{}/{} files; {:.1f} total mins elapsed' + .format(ifn, nfiles, (time()-start)/60.)) + urat = fitsio.read(fn) + idurat, idobjs, dist = radec_match_to(urat, objs, sep=matchrad, + return_sep=True) + + # ADM update matches whenever we have a CLOSER match. + ii = (distance[idobjs] == -1) | (distance[idobjs] > dist) + done[idobjs[ii]] = urat[idurat[ii]] + distance[idobjs[ii]] = dist[ii] + + # ADM add the distances to the output array. + dt = uratdatamodel.dtype.descr + [("DISTANCE", ">f4")] + output = np.zeros(nobjs, dtype=dt) + for col in uratdatamodel.dtype.names: + output[col] = done[col] + output["DISTANCE"] = distance + + return output From 78517bcf3302d4235f52d596e9d18b1fd70a7c9e Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Wed, 31 Jul 2019 16:22:33 -0700 Subject: [PATCH 11/30] update magnitude limit for GFAs to 21 --- bin/select_gfas | 4 ++-- py/desitarget/gfa.py | 6 +++--- py/desitarget/uratmatch.py | 16 ++++++++-------- 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/bin/select_gfas b/bin/select_gfas index 9067bc81f..95c87c834 100755 --- a/bin/select_gfas +++ b/bin/select_gfas @@ -23,8 +23,8 @@ ap.add_argument("surveydir", ap.add_argument("dest", help="Output GFA targets file (e.g. '/project/projectdirs/desi/target/catalogs/gfas-dr6-0.20.0.fits' at NERSC)") ap.add_argument('-m', '--maglim', type=float, - help="Magnitude limit on GFA targets in Gaia G-band (defaults to [18])", - default=18) + help="Magnitude limit on GFA targets in Gaia G-band (defaults to [21])", + default=21) ap.add_argument('-n', "--numproc", type=int, help="number of concurrent processes to use (defaults to [{}])".format(nproc), default=nproc) diff --git a/py/desitarget/gfa.py b/py/desitarget/gfa.py index e86b686cb..35a006fc1 100644 --- a/py/desitarget/gfa.py +++ b/py/desitarget/gfa.py @@ -404,9 +404,9 @@ def _update_status(result): # ADM a final clean-up to remove columns that are NaN (from # ADM Gaia-matching) or that are exactly 0 (in the sweeps). - for col in ["PMRA", "PMDEC"]: - ii = ~np.isnan(gfas[col]) & (gfas[col] != 0) - gfas = gfas[ii] +# for col in ["PMRA", "PMDEC"]: +# ii = ~np.isnan(gfas[col]) & (gfas[col] != 0) +# gfas = gfas[ii] # ADM limit to DESI footprint or passed tiles, if not cmx'ing. if not cmx: diff --git a/py/desitarget/uratmatch.py b/py/desitarget/uratmatch.py index 90f273776..c4501cc6d 100644 --- a/py/desitarget/uratmatch.py +++ b/py/desitarget/uratmatch.py @@ -179,7 +179,7 @@ def urat_csv_to_fits(numproc=5): ---------- numproc : :class:`int`, optional, defaults to 5 The number of parallel processes to use. - + Returns ------- Nothing @@ -197,7 +197,7 @@ def urat_csv_to_fits(numproc=5): """ # ADM the resolution at which the URAT HEALPix files should be stored. nside = _get_urat_nside() - + # ADM check that the URAT_DIR is set. uratdir = _get_urat_dir() log.info("running on {} processors".format(numproc)) @@ -469,7 +469,7 @@ def make_urat_files(numproc=5, download=False): urat_binary_to_csv(numproc=numproc) log.info('Converted binary files to CSV...t={:.1f}s'.format(time()-t0)) - + urat_csv_to_fits(numproc=numproc) log.info('Converted CSV files to FITS...t={:.1f}s'.format(time()-t0)) @@ -533,7 +533,7 @@ def find_urat_files(objs, neighbors=True, radec=False): # ADM restrict to only files/HEALPixels actually covered by URAT. uratfiles = [fn for fn in uratfiles if os.path.exists(fn)] - + return uratfiles @@ -570,11 +570,11 @@ def match_to_urat(objs, matchrad=1.): # ADM objects without matches should have URAT_ID, DISTANCE of -1. done["URAT_ID"] = -1 distance = np.zeros(nobjs) - 1 - + # ADM determine which URAT files need to be scraped. uratfiles = find_urat_files(objs) nfiles = len(uratfiles) - + # ADM catch the case of no matches to URAT. if nfiles == 0: return done @@ -592,12 +592,12 @@ def match_to_urat(objs, matchrad=1.): ii = (distance[idobjs] == -1) | (distance[idobjs] > dist) done[idobjs[ii]] = urat[idurat[ii]] distance[idobjs[ii]] = dist[ii] - + # ADM add the distances to the output array. dt = uratdatamodel.dtype.descr + [("DISTANCE", ">f4")] output = np.zeros(nobjs, dtype=dt) for col in uratdatamodel.dtype.names: output[col] = done[col] output["DISTANCE"] = distance - + return output From eb8ebc550f15f01efc2c7b781ed7258e04889392 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Wed, 31 Jul 2019 16:38:44 -0700 Subject: [PATCH 12/30] add extra keywords to the header of the GFAs file --- bin/select_gfas | 7 ++++++- py/desitarget/io.py | 18 +++++++++++++----- 2 files changed, 19 insertions(+), 6 deletions(-) diff --git a/bin/select_gfas b/bin/select_gfas index 95c87c834..76778508e 100755 --- a/bin/select_gfas +++ b/bin/select_gfas @@ -76,6 +76,11 @@ log.info('running on {} processors...t = {:.1f} mins'.format(ns.numproc, (time() gfas = select_gfas(infiles, maglim=ns.maglim, numproc=ns.numproc, tilesfile=ns.tiles, cmx=cmx, mindec=ns.mindec, mingalb=ns.mingalb) -io.write_gfas(ns.dest, gfas, indir=ns.surveydir, indir2=ns.surveydir2, nside=nside, survey=survey) +# ADM extra header keywords for the output fits file. +extra = {k: v for k, v in zip(["maglim", "mindec", "mingalb"], + [ns.maglim, ns.mindec, ns.mingalb])} + +io.write_gfas(ns.dest, gfas, indir=ns.surveydir, indir2=ns.surveydir2, + nside=nside, survey=survey, extra=extra) log.info('{} GFAs written to {}...t = {:.1f} mins'.format(len(gfas), ns.dest, (time()-t0)/60.)) diff --git a/py/desitarget/io.py b/py/desitarget/io.py index b1b632b0c..9a5880d70 100644 --- a/py/desitarget/io.py +++ b/py/desitarget/io.py @@ -554,7 +554,7 @@ def write_skies(filename, data, indir=None, indir2=None, supp=False, def write_gfas(filename, data, indir=None, indir2=None, nside=None, - survey="?"): + survey="?", extra=None): """Write a catalogue of Guide/Focus/Alignment targets. Parameters @@ -564,13 +564,16 @@ def write_gfas(filename, data, indir=None, indir2=None, nside=None, data : :class:`~numpy.ndarray` Array of GFAs to write to file. indir, indir2 : :class:`str`, optional, defaults to None. - Name of input Legacy Survey Data Release directory or directories, - write to header of output file if passed (and if not None). + Legacy Survey Data Release directory or directories, write to + header of output file if passed (and if not None). nside: :class:`int`, defaults to None. - If passed, add a column to the GFAs array popluated with HEALPixels - at resolution `nside`. + If passed, add a column to the GFAs array popluated with + HEALPixels at resolution `nside`. survey : :class:`str`, optional, defaults to "?" Written to output file header as the keyword `SURVEY`. + extra :class:`dict`, optional + If passed (and not None), write these extra dictionary keys and + values to the output header. """ # ADM rename 'TYPE' to 'MORPHTYPE'. data = rfn.rename_fields(data, {'TYPE': 'MORPHTYPE'}) @@ -590,6 +593,11 @@ def write_gfas(filename, data, indir=None, indir2=None, nside=None, if indir2 is not None: depend.setdep(hdr, 'input-data-release-2', indir2) + # ADM add the extra dictionary to the header. + if extra is not None: + for key in extra: + hdr[key] = extra[key] + # ADM add HEALPix column, if requested by input. if nside is not None: theta, phi = np.radians(90-data["DEC"]), np.radians(data["RA"]) From d87c521d67d85e429caef6f559c02d14137d81b8 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Wed, 31 Jul 2019 17:24:05 -0700 Subject: [PATCH 13/30] enforce Dec range in hp_in_box, even when returning the full RA range --- py/desitarget/geomask.py | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/py/desitarget/geomask.py b/py/desitarget/geomask.py index 74312ee78..286972e04 100644 --- a/py/desitarget/geomask.py +++ b/py/desitarget/geomask.py @@ -853,33 +853,33 @@ def hp_in_box(nside, radecbox, inclusive=True, fact=4): - When the RA range exceeds 180o, `healpy.query_polygon()` defines the range as that with the smallest area (i.e the box can wrap-around in RA). To avoid any ambiguity, this function - will return ALL HEALPixels in such cases. + will only limit by the passed Decs in such cases. - Only strictly correct for Decs from -90+1e-5(o) to 90-1e5(o). """ ramin, ramax, decmin, decmax = radecbox # ADM area enclosed isn't well-defined if RA covers more than 180o. - if np.abs(ramax-ramin) > 180.: + if np.abs(ramax-ramin) <= 180.: + # ADM retrieve RA range. The 1e-5 prevents edge effects near poles. + npole, spole = 90-1e-5, -90+1e-5 + # ADM convert RA/Dec to co-latitude and longitude in radians. + rapairs = np.array([ramin, ramin, ramax, ramax]) + decpairs = np.array([spole, npole, npole, spole]) + thetapairs, phipairs = np.radians(90.-decpairs), np.radians(rapairs) + + # ADM convert to Cartesian vectors remembering to transpose + # ADM to pass the array to query_polygon in the correct order. + vecs = hp.dir2vec(thetapairs, phipairs).T + + # ADM determine the pixels that touch the RA range. + pixra = hp.query_polygon(nside, vecs, + inclusive=inclusive, fact=fact, nest=True) + else: log.warning('Max RA ({}) and Min RA ({}) separated by > 180o...' .format(ramax, ramin)) - log.warning('...returning full set of HEALPixels at nside={}' + log.warning('...will only limit to passed Declinations' .format(nside)) - return np.arange(hp.nside2npix(nside)) - - # ADM retrieve RA range. The 1e-5 prevents edge effects near poles. - npole, spole = 90-1e-5, -90+1e-5 - # ADM convert RA/Dec to co-latitude and longitude in radians. - rapairs = np.array([ramin, ramin, ramax, ramax]) - decpairs = np.array([spole, npole, npole, spole]) - thetapairs, phipairs = np.radians(90.-decpairs), np.radians(rapairs) - - # ADM convert to Cartesian vectors remembering to transpose - # ADM to pass the array to query_polygon in the correct order. - vecs = hp.dir2vec(thetapairs, phipairs).T - - # ADM determine the pixels that touch the RA range. - pixra = hp.query_polygon(nside, vecs, - inclusive=inclusive, fact=fact, nest=True) + pixra = np.arange(hp.nside2npix(nside)) # ADM determine the pixels that touch the Dec range. pixdec = hp_in_dec_range(nside, decmin, decmax, inclusive=inclusive) @@ -928,7 +928,7 @@ def hp_in_dec_range(nside, decmin, decmax, inclusive=True): def hp_beyond_gal_b(nside, mingalb, neighbors=True): - """Find all HEALPixels with centers and neigbors beyond a Galactic b. + """Find HEALPixels with centers and neighbors beyond a Galactic b. Parameters ---------- From 80c46f94c06b09abb41132f00fd4beec54f70ffc Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Wed, 31 Jul 2019 18:01:43 -0700 Subject: [PATCH 14/30] better limit the Gaia files being considered when making GFAs --- py/desitarget/gfa.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/py/desitarget/gfa.py b/py/desitarget/gfa.py index 35a006fc1..f52794faa 100644 --- a/py/desitarget/gfa.py +++ b/py/desitarget/gfa.py @@ -17,7 +17,7 @@ import desitarget.io from desitarget.internal import sharedmem -from desitarget.gaiamatch import read_gaia_file +from desitarget.gaiamatch import read_gaia_file, find_gaia_files_beyond_gal_b from desitarget.gaiamatch import find_gaia_files_tiles, find_gaia_files_box from desitarget.targets import encode_targetid, resolve from desitarget.geomask import is_in_gal_box, is_in_box @@ -251,7 +251,9 @@ def all_gaia_in_tiles(maglim=18, numproc=4, allsky=False, """ # ADM grab paths to Gaia files in the sky or the DESI footprint. if allsky: - infiles = find_gaia_files_box([0, 360, mindec, 90]) + infilesbox = find_gaia_files_box([0, 360, mindec, 90]) + infilesgalb = find_gaia_files_beyond_gal_b(mingalb) + infiles = list(set(infilesbox).intersection(set(infilesgalb))) else: infiles = find_gaia_files_tiles(tiles=tiles, neighbors=False) nfiles = len(infiles) From 97252a4c0664bcaec3504006bc8170adeb8ebbd0 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 10:29:27 -0700 Subject: [PATCH 15/30] parallelized function to match URAT to zero-PM GFAs; limit to numproc=4 where appropriate --- bin/select_gfas | 1 + py/desitarget/gfa.py | 122 ++++++++++++++++++++++++++++++++++--- py/desitarget/uratmatch.py | 46 ++++++++------ 3 files changed, 140 insertions(+), 29 deletions(-) diff --git a/bin/select_gfas b/bin/select_gfas index 76778508e..9d3ede485 100755 --- a/bin/select_gfas +++ b/bin/select_gfas @@ -80,6 +80,7 @@ gfas = select_gfas(infiles, maglim=ns.maglim, numproc=ns.numproc, tilesfile=ns.t extra = {k: v for k, v in zip(["maglim", "mindec", "mingalb"], [ns.maglim, ns.mindec, ns.mingalb])} +log.info('Writing GFAs to file...t = {:.1f} mins'.format(len(gfas), ns.dest, (time()-t0)/60.)) io.write_gfas(ns.dest, gfas, indir=ns.surveydir, indir2=ns.surveydir2, nside=nside, survey=survey, extra=extra) diff --git a/py/desitarget/gfa.py b/py/desitarget/gfa.py index f52794faa..a0ef5da3a 100644 --- a/py/desitarget/gfa.py +++ b/py/desitarget/gfa.py @@ -29,7 +29,6 @@ bricks = brick.Bricks(bricksize=0.25) # ADM set up the default DESI logger. log = get_logger() -start = time() # ADM the current data model for columns in the GFA files. gfadatamodel = np.array([], dtype=[ @@ -45,7 +44,7 @@ ('GAIA_PHOT_G_MEAN_MAG', '>f4'), ('GAIA_PHOT_G_MEAN_FLUX_OVER_ERROR', '>f4'), ('GAIA_PHOT_BP_MEAN_MAG', '>f4'), ('GAIA_PHOT_BP_MEAN_FLUX_OVER_ERROR', '>f4'), ('GAIA_PHOT_RP_MEAN_MAG', '>f4'), ('GAIA_PHOT_RP_MEAN_FLUX_OVER_ERROR', '>f4'), - ('GAIA_ASTROMETRIC_EXCESS_NOISE', '>f4') + ('GAIA_ASTROMETRIC_EXCESS_NOISE', '>f4'), ('URAT_ID', '>i8') ]) @@ -127,7 +126,7 @@ def gaia_gfas_from_sweep(filename, maglim=18.): # ADM remove the TARGETID, BRICK_OBJID, REF_CAT, REF_EPOCH columns # ADM and populate them later as they require special treatment. cols = list(gfadatamodel.dtype.names) - for col in ["TARGETID", "BRICK_OBJID", "REF_CAT", "REF_EPOCH"]: + for col in ["TARGETID", "BRICK_OBJID", "REF_CAT", "REF_EPOCH", "URAT_ID"]: cols.remove(col) for col in cols: gfas[col] = objects[col] @@ -297,6 +296,92 @@ def _update_status(result): return gfas +def add_urat_pms(objs, numproc=4): + """Add proper motions from URAT to a set of objects. + + Parameters + ---------- + objs : :class:`~numpy.ndarray` + Array of objects to update. Must include the columns "PMRA", + "PMDEC", "REF_ID" (unique for each object) and "URAT_ID". + numproc : :class:`int`, optional, defaults to 4 + The number of parallel processes to use. + + Returns + ------- + :class:`~numpy.ndarray` + The input array with the "PMRA", PMDEC" and "URAT_ID" + columns updated to include URAT information. + + Notes + ----- + - Order is retained using "REF_ID": The input and output + arrays should have the same order. + """ + # ADM check REF_ID is indeed unique for each object. + assert len(objs["REF_ID"]) == len(np.unique(objs["REF_ID"])) + + # ADM record the original REF_IDs so we can match back to them. + origids = objs["REF_ID"] + + # ADM loosely group the input objects on the sky. The + # ADM choice of NSIDE=2 isn't magical, the URAT-matching + # ADM code is just quicker fo clumped objects. + theta, phi = np.radians(90-objs["DEC"]), np.radians(objs["RA"]) + pixels = hp.ang2pix(2, theta, phi, nest=True) + + # ADM reorder objects (and pixels themselves) based on pixel number. + ii = np.argsort(pixels) + objs, pixels = objs[ii], pixels[ii] + + # ADM create pixel-split sub-lists of the objects. + # ADM here, np.diff marks the transition to the next pixel number. + splitobjs = np.split(objs, np.where(np.diff(pixels))[0]+1) + nallpix = len(splitobjs) + + # ADM function to run on each of the HEALPix-split input objs. + def _get_urat_matches(splitobj): + '''wrapper on match_to_urat() for rec array (matchrad=0.1")''' + # ADM also return the REF_ID to track the objects. + return [match_to_urat(splitobj, matchrad=0.1), splitobj["REF_ID"]] + + # ADM this is just to count pixels in _update_status. + npix = np.zeros((), dtype='i8') + t0 = time() + + def _update_status(result): + """wrapper function for the critical reduction operation, + that occurs on the main parallel process""" + if npix % 5 == 0 and npix > 0: + elapsed = (time()-t0)/60. + rate = 60*elapsed/npix + log.info('{}/{} pixels; {:.1f} sec/pix...t = {:.1f} mins' + .format(npix, nallpix, rate, elapsed)) + npix[...] += 1 # this is an in-place modification. + return result + + # - Parallel process pixels. + if numproc > 1: + pool = sharedmem.MapReduce(np=numproc) + with pool: + urats = pool.map(_get_urat_matches, splitobjs, reduce=_update_status) + else: + urats = [] + for splitobj in splitobjs: + urats.append(_update_status(_get_urat_matches(splitobj))) + + # ADM remember to grab the REFIDs as well as the URAT matches. + refids = np.concatenate(np.array(urats)[:,1]) + urats = np.concatenate(np.array(urats)[:,0]) + + # ADM sort the output to match the input, on REF_ID. + ii = np.zeros_like(refids) + ii[np.argsort(origids)] = np.argsort(refids) + assert np.all(refids[ii] == origids) + + return urats[ii] + + def select_gfas(infiles, maglim=18, numproc=4, tilesfile=None, cmx=False, mindec=-30, mingalb=10): """Create a set of GFA locations using Gaia and matching to sweeps. @@ -336,6 +421,12 @@ def select_gfas(infiles, maglim=18, numproc=4, tilesfile=None, - The tiles loaded from `tilesfile` will only be those in DESI. So, for custom tilings, set IN_DESI==1 in your tiles file. """ + # ADM force to no more than numproc=4 for I/O heavy processes. + numproc4 = numproc + if numproc4 > 4: + log.info('Forcing numproc to 4 for I/O heavy parts of code') + numproc4 = 4 + # ADM convert a single file, if passed to a list of files. if isinstance(infiles, str): infiles = [infiles, ] @@ -362,8 +453,8 @@ def _get_gfas(fn): return gaia_gfas_from_sweep(fn, maglim=maglim) # ADM this is just to count sweeps files in _update_status. - nfile = np.zeros((), dtype='i8') t0 = time() + nfile = np.zeros((), dtype='i8') def _update_status(result): """wrapper function for the critical reduction operation, @@ -377,8 +468,8 @@ def _update_status(result): return result # - Parallel process input files. - if numproc > 1: - pool = sharedmem.MapReduce(np=numproc) + if numproc4 > 1: + pool = sharedmem.MapReduce(np=numproc4) with pool: gfas = pool.map(_get_gfas, infiles, reduce=_update_status) else: @@ -394,7 +485,7 @@ def _update_status(result): # ADM retrieve Gaia objects in the DESI footprint or passed tiles. log.info('Retrieving additional Gaia objects...t = {:.1f} mins' .format((time()-t0)/60)) - gaia = all_gaia_in_tiles(maglim=maglim, numproc=numproc, allsky=cmx, + gaia = all_gaia_in_tiles(maglim=maglim, numproc=numproc4, allsky=cmx, tiles=tiles, mindec=mindec, mingalb=mingalb) # ADM remove any duplicates. Order is important here, as np.unique @@ -404,11 +495,22 @@ def _update_status(result): _, ind = np.unique(gfas["REF_ID"], return_index=True) gfas = gfas[ind] + # ADM for zero/NaN proper motion objects, add in URAT proper motions. + ii = ((np.isnan(gfas["PMRA"]) | (gfas["PMRA"] == 0)) & + (np.isnan(gfas["PMDEC"]) | (gfas["PMDEC"] == 0))) + log.info('Adding URAT for {} objects with no proper motion...t = {:.1f} mins' + .format(np.sum(ii), (time()-t0)/60)) + urat = add_urat_pms(gfas[ii], numproc=numproc) + log.info('Found an additional {} URAT objects...t = {:.1f} mins' + .format(np.sum(urat["REF_ID"] != -1), (time()-t0)/60)) + for col in "PMRA", "PMDEC", "URAT_ID": + gfas[col][ii] = urat[col] + # ADM a final clean-up to remove columns that are NaN (from # ADM Gaia-matching) or that are exactly 0 (in the sweeps). -# for col in ["PMRA", "PMDEC"]: -# ii = ~np.isnan(gfas[col]) & (gfas[col] != 0) -# gfas = gfas[ii] + ii = ((np.isnan(gfas["PMRA"]) | (gfas["PMRA"] == 0)) & + (np.isnan(gfas["PMDEC"]) | (gfas["PMDEC"] == 0))) + gfas = gfas[ii] # ADM limit to DESI footprint or passed tiles, if not cmx'ing. if not cmx: diff --git a/py/desitarget/uratmatch.py b/py/desitarget/uratmatch.py index c4501cc6d..e3af55509 100644 --- a/py/desitarget/uratmatch.py +++ b/py/desitarget/uratmatch.py @@ -537,7 +537,7 @@ def find_urat_files(objs, neighbors=True, radec=False): return uratfiles -def match_to_urat(objs, matchrad=1.): +def match_to_urat(objs, matchrad=1., radec=False): """Match objects to URAT healpix files and return URAT information. Parameters @@ -546,6 +546,9 @@ def match_to_urat(objs, matchrad=1.): Must contain at least "RA" and "DEC". matchrad : :class:`float`, optional, defaults to 1 arcsec The radius at which to match in arcseconds. + radec : :class:`bool`, optional, defaults to ``False`` + If ``True`` then the passed `objs` is an [RA, Dec] list instead of + a rec array. Returns ------- @@ -563,8 +566,14 @@ def match_to_urat(objs, matchrad=1.): - Because this reads in HEALPixel split files, it's (far) faster for objects that are clumped rather than widely distributed. """ + # ADM parse whether a structure or coordinate list was passed. + if radec: + ra, dec = objs + else: + ra, dec = objs["RA"], objs["DEC"] + # ADM set up an array of URAT information for the output. - nobjs = len(objs) + nobjs = len(ra) done = np.zeros(nobjs, dtype=uratdatamodel.dtype) # ADM objects without matches should have URAT_ID, DISTANCE of -1. @@ -572,26 +581,25 @@ def match_to_urat(objs, matchrad=1.): distance = np.zeros(nobjs) - 1 # ADM determine which URAT files need to be scraped. - uratfiles = find_urat_files(objs) + uratfiles = find_urat_files([ra, dec], radec=True) nfiles = len(uratfiles) # ADM catch the case of no matches to URAT. - if nfiles == 0: - return done - - # ADM loop through the URAT files and find matches. - for ifn, fn in enumerate(uratfiles): - if ifn % 100 == 0 and ifn > 0: - log.info('{}/{} files; {:.1f} total mins elapsed' - .format(ifn, nfiles, (time()-start)/60.)) - urat = fitsio.read(fn) - idurat, idobjs, dist = radec_match_to(urat, objs, sep=matchrad, - return_sep=True) - - # ADM update matches whenever we have a CLOSER match. - ii = (distance[idobjs] == -1) | (distance[idobjs] > dist) - done[idobjs[ii]] = urat[idurat[ii]] - distance[idobjs[ii]] = dist[ii] + if nfiles > 0: + # ADM loop through the URAT files and find matches. + for ifn, fn in enumerate(uratfiles): + if ifn % 500 == 0 and ifn > 0: + log.info('{}/{} files; {:.1f} total mins elapsed' + .format(ifn, nfiles, (time()-start)/60.)) + urat = fitsio.read(fn) + idurat, idobjs, dist = radec_match_to( + [urat["RA"], urat["DEC"]], [ra, dec], + sep=matchrad, radec=True, return_sep=True) + + # ADM update matches whenever we have a CLOSER match. + ii = (distance[idobjs] == -1) | (distance[idobjs] > dist) + done[idobjs[ii]] = urat[idurat[ii]] + distance[idobjs[ii]] = dist[ii] # ADM add the distances to the output array. dt = uratdatamodel.dtype.descr + [("DISTANCE", ">f4")] From d9aeb79719f88a6167840595bb9fa937ca33680c Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 10:48:18 -0700 Subject: [PATCH 16/30] code style and extra imports --- py/desitarget/gfa.py | 6 ++++-- py/desitarget/randoms.py | 2 +- py/desitarget/uratmatch.py | 5 +++-- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/py/desitarget/gfa.py b/py/desitarget/gfa.py index a0ef5da3a..30f6850e3 100644 --- a/py/desitarget/gfa.py +++ b/py/desitarget/gfa.py @@ -10,6 +10,7 @@ import glob import os from time import time +import healpy as hp import desimodel.focalplane import desimodel.io @@ -19,6 +20,7 @@ from desitarget.internal import sharedmem from desitarget.gaiamatch import read_gaia_file, find_gaia_files_beyond_gal_b from desitarget.gaiamatch import find_gaia_files_tiles, find_gaia_files_box +from desitarget.uratmatch import match_to_urat from desitarget.targets import encode_targetid, resolve from desitarget.geomask import is_in_gal_box, is_in_box @@ -371,8 +373,8 @@ def _update_status(result): urats.append(_update_status(_get_urat_matches(splitobj))) # ADM remember to grab the REFIDs as well as the URAT matches. - refids = np.concatenate(np.array(urats)[:,1]) - urats = np.concatenate(np.array(urats)[:,0]) + refids = np.concatenate(np.array(urats)[:, 1]) + urats = np.concatenate(np.array(urats)[:, 0]) # ADM sort the output to match the input, on REF_ID. ii = np.zeros_like(refids) diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py index faa9b7e05..6b9d0d98b 100644 --- a/py/desitarget/randoms.py +++ b/py/desitarget/randoms.py @@ -751,7 +751,7 @@ def stellar_density(nside=256): pixout = np.zeros(npix, dtype='int32') # ADM find all of the Gaia files. - filenames = glob(os.path.join(hpdir,'*fits')) + filenames = glob(os.path.join(hpdir, '*fits')) # ADM read in each file, restricting to the criteria for point # ADM sources and storing in a HEALPixel map at resolution nside. diff --git a/py/desitarget/uratmatch.py b/py/desitarget/uratmatch.py index e3af55509..7b2787784 100644 --- a/py/desitarget/uratmatch.py +++ b/py/desitarget/uratmatch.py @@ -41,6 +41,7 @@ ('PMRA', '>f4'), ('PMDEC', '>f4'), ('PM_ERROR', '>f4') ]) + def _get_urat_dir(): """Convenience function to grab the URAT environment variable. @@ -560,7 +561,7 @@ def match_to_urat(objs, matchrad=1., radec=False): Notes ----- - - For objects that do NOT have a match in URAT, the "URAT_ID" + - For objects that do NOT have a match in URAT, the "URAT_ID" and "DISTANCE" columns are -1, and other columns are zero. - Retrieves the CLOSEST match to URAT for each passed object. - Because this reads in HEALPixel split files, it's (far) faster @@ -590,7 +591,7 @@ def match_to_urat(objs, matchrad=1., radec=False): for ifn, fn in enumerate(uratfiles): if ifn % 500 == 0 and ifn > 0: log.info('{}/{} files; {:.1f} total mins elapsed' - .format(ifn, nfiles, (time()-start)/60.)) + .format(ifn, nfiles, (time()-start)/60.)) urat = fitsio.read(fn) idurat, idobjs, dist = radec_match_to( [urat["RA"], urat["DEC"]], [ra, dec], From 93730116f583cece4e257820deb635b8c4db6abf Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 12:04:32 -0700 Subject: [PATCH 17/30] optimize parallelization of matching; add URAT matching distance to output array --- bin/select_gfas | 2 +- py/desitarget/gfa.py | 39 ++++++++++++++++++++------------------ py/desitarget/uratmatch.py | 18 +++++++++--------- 3 files changed, 31 insertions(+), 28 deletions(-) diff --git a/bin/select_gfas b/bin/select_gfas index 9d3ede485..c685a429f 100755 --- a/bin/select_gfas +++ b/bin/select_gfas @@ -26,7 +26,7 @@ ap.add_argument('-m', '--maglim', type=float, help="Magnitude limit on GFA targets in Gaia G-band (defaults to [21])", default=21) ap.add_argument('-n', "--numproc", type=int, - help="number of concurrent processes to use (defaults to [{}])".format(nproc), + help="number of concurrent processes to use (defaults to [{}]). Note that if numproc > 4, I/O limited parts of the code revert to numproc=4".format(nproc), default=nproc) ap.add_argument('-s2', "--surveydir2", help='Additional Legacy Surveys Data Release directory (useful for combining, e.g., DR8 into one file of GFAs)', diff --git a/py/desitarget/gfa.py b/py/desitarget/gfa.py index 30f6850e3..0d47adb74 100644 --- a/py/desitarget/gfa.py +++ b/py/desitarget/gfa.py @@ -46,7 +46,7 @@ ('GAIA_PHOT_G_MEAN_MAG', '>f4'), ('GAIA_PHOT_G_MEAN_FLUX_OVER_ERROR', '>f4'), ('GAIA_PHOT_BP_MEAN_MAG', '>f4'), ('GAIA_PHOT_BP_MEAN_FLUX_OVER_ERROR', '>f4'), ('GAIA_PHOT_RP_MEAN_MAG', '>f4'), ('GAIA_PHOT_RP_MEAN_FLUX_OVER_ERROR', '>f4'), - ('GAIA_ASTROMETRIC_EXCESS_NOISE', '>f4'), ('URAT_ID', '>i8') + ('GAIA_ASTROMETRIC_EXCESS_NOISE', '>f4'), ('URAT_ID', '>i8'), ('URAT_SEP', '>f4') ]) @@ -128,7 +128,8 @@ def gaia_gfas_from_sweep(filename, maglim=18.): # ADM remove the TARGETID, BRICK_OBJID, REF_CAT, REF_EPOCH columns # ADM and populate them later as they require special treatment. cols = list(gfadatamodel.dtype.names) - for col in ["TARGETID", "BRICK_OBJID", "REF_CAT", "REF_EPOCH", "URAT_ID"]: + for col in ["TARGETID", "BRICK_OBJID", "REF_CAT", "REF_EPOCH", + "URAT_ID", "URAT_SEP"]: cols.remove(col) for col in cols: gfas[col] = objects[col] @@ -305,14 +306,14 @@ def add_urat_pms(objs, numproc=4): ---------- objs : :class:`~numpy.ndarray` Array of objects to update. Must include the columns "PMRA", - "PMDEC", "REF_ID" (unique for each object) and "URAT_ID". + "PMDEC", "REF_ID" (unique per object) "URAT_ID" and "URAT_SEP". numproc : :class:`int`, optional, defaults to 4 The number of parallel processes to use. Returns ------- :class:`~numpy.ndarray` - The input array with the "PMRA", PMDEC" and "URAT_ID" + The input array with the "PMRA", PMDEC", "URAT_ID" and "URAT_SEP" columns updated to include URAT information. Notes @@ -326,11 +327,11 @@ def add_urat_pms(objs, numproc=4): # ADM record the original REF_IDs so we can match back to them. origids = objs["REF_ID"] - # ADM loosely group the input objects on the sky. The - # ADM choice of NSIDE=2 isn't magical, the URAT-matching - # ADM code is just quicker fo clumped objects. + # ADM loosely group the input objects on the sky. NSIDE=16 seems + # ADM to nicely balance sample sizes for matching, with the code + # ADM being quicker for clumped objects because of file I/O. theta, phi = np.radians(90-objs["DEC"]), np.radians(objs["RA"]) - pixels = hp.ang2pix(2, theta, phi, nest=True) + pixels = hp.ang2pix(16, theta, phi, nest=True) # ADM reorder objects (and pixels themselves) based on pixel number. ii = np.argsort(pixels) @@ -343,9 +344,9 @@ def add_urat_pms(objs, numproc=4): # ADM function to run on each of the HEALPix-split input objs. def _get_urat_matches(splitobj): - '''wrapper on match_to_urat() for rec array (matchrad=0.1")''' + '''wrapper on match_to_urat() for rec array (matchrad=0.5")''' # ADM also return the REF_ID to track the objects. - return [match_to_urat(splitobj, matchrad=0.1), splitobj["REF_ID"]] + return [match_to_urat(splitobj, matchrad=0.5), splitobj["REF_ID"]] # ADM this is just to count pixels in _update_status. npix = np.zeros((), dtype='i8') @@ -354,10 +355,10 @@ def _get_urat_matches(splitobj): def _update_status(result): """wrapper function for the critical reduction operation, that occurs on the main parallel process""" - if npix % 5 == 0 and npix > 0: + if npix % 200 == 0 and npix > 0: elapsed = (time()-t0)/60. - rate = 60*elapsed/npix - log.info('{}/{} pixels; {:.1f} sec/pix...t = {:.1f} mins' + rate = npix/elapsed/60. + log.info('{}/{} pixels; {:.1f} pix/sec...t = {:.1f} mins' .format(npix, nallpix, rate, elapsed)) npix[...] += 1 # this is an in-place modification. return result @@ -419,14 +420,16 @@ def select_gfas(infiles, maglim=18, numproc=4, tilesfile=None, Notes ----- - - If numproc==1, use the serial code instead of the parallel code. + - If numproc==1, use the serial code instead of parallel code. + - If numproc > 4, then numproc=4 is enforced for (just those) + parts of the code that are I/O limited. - The tiles loaded from `tilesfile` will only be those in DESI. So, for custom tilings, set IN_DESI==1 in your tiles file. """ - # ADM force to no more than numproc=4 for I/O heavy processes. + # ADM force to no more than numproc=4 for I/O limited processes. numproc4 = numproc if numproc4 > 4: - log.info('Forcing numproc to 4 for I/O heavy parts of code') + log.info('Forcing numproc to 4 for I/O limited parts of code') numproc4 = 4 # ADM convert a single file, if passed to a list of files. @@ -504,8 +507,8 @@ def _update_status(result): .format(np.sum(ii), (time()-t0)/60)) urat = add_urat_pms(gfas[ii], numproc=numproc) log.info('Found an additional {} URAT objects...t = {:.1f} mins' - .format(np.sum(urat["REF_ID"] != -1), (time()-t0)/60)) - for col in "PMRA", "PMDEC", "URAT_ID": + .format(np.sum(urat["URAT_ID"] != -1), (time()-t0)/60)) + for col in "PMRA", "PMDEC", "URAT_ID", "URAT_SEP": gfas[col][ii] = urat[col] # ADM a final clean-up to remove columns that are NaN (from diff --git a/py/desitarget/uratmatch.py b/py/desitarget/uratmatch.py index 7b2787784..8f383352f 100644 --- a/py/desitarget/uratmatch.py +++ b/py/desitarget/uratmatch.py @@ -556,13 +556,13 @@ def match_to_urat(objs, matchrad=1., radec=False): :class:`~numpy.ndarray` The matching URAT information for each object. The returned format is as for desitarget.uratmatch.uratdatamodel with - and extra column "DISTANCE" which is the matching distance + and extra column "URAT_SEP" which is the matching distance in ARCSECONDS. Notes ----- - For objects that do NOT have a match in URAT, the "URAT_ID" - and "DISTANCE" columns are -1, and other columns are zero. + and "URAT_SEP" columns are -1, and other columns are zero. - Retrieves the CLOSEST match to URAT for each passed object. - Because this reads in HEALPixel split files, it's (far) faster for objects that are clumped rather than widely distributed. @@ -577,9 +577,9 @@ def match_to_urat(objs, matchrad=1., radec=False): nobjs = len(ra) done = np.zeros(nobjs, dtype=uratdatamodel.dtype) - # ADM objects without matches should have URAT_ID, DISTANCE of -1. + # ADM objects without matches should have URAT_ID, URAT_SEP of -1. done["URAT_ID"] = -1 - distance = np.zeros(nobjs) - 1 + urat_sep = np.zeros(nobjs) - 1 # ADM determine which URAT files need to be scraped. uratfiles = find_urat_files([ra, dec], radec=True) @@ -598,15 +598,15 @@ def match_to_urat(objs, matchrad=1., radec=False): sep=matchrad, radec=True, return_sep=True) # ADM update matches whenever we have a CLOSER match. - ii = (distance[idobjs] == -1) | (distance[idobjs] > dist) + ii = (urat_sep[idobjs] == -1) | (urat_sep[idobjs] > dist) done[idobjs[ii]] = urat[idurat[ii]] - distance[idobjs[ii]] = dist[ii] + urat_sep[idobjs[ii]] = dist[ii] - # ADM add the distances to the output array. - dt = uratdatamodel.dtype.descr + [("DISTANCE", ">f4")] + # ADM add the separation distances to the output array. + dt = uratdatamodel.dtype.descr + [("URAT_SEP", ">f4")] output = np.zeros(nobjs, dtype=dt) for col in uratdatamodel.dtype.names: output[col] = done[col] - output["DISTANCE"] = distance + output["URAT_SEP"] = urat_sep return output From 0567210e64953f51ef880b89d99d4cf1472dcb8c Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 12:11:57 -0700 Subject: [PATCH 18/30] debugging --- py/desitarget/gfa.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/py/desitarget/gfa.py b/py/desitarget/gfa.py index 0d47adb74..dc625c430 100644 --- a/py/desitarget/gfa.py +++ b/py/desitarget/gfa.py @@ -515,7 +515,7 @@ def _update_status(result): # ADM Gaia-matching) or that are exactly 0 (in the sweeps). ii = ((np.isnan(gfas["PMRA"]) | (gfas["PMRA"] == 0)) & (np.isnan(gfas["PMDEC"]) | (gfas["PMDEC"] == 0))) - gfas = gfas[ii] + gfas = gfas[~ii] # ADM limit to DESI footprint or passed tiles, if not cmx'ing. if not cmx: From c8f17c29e99dc8bdd26f0caf901280db5b91d84a Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 12:20:34 -0700 Subject: [PATCH 19/30] update changes docs --- doc/changes.rst | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/doc/changes.rst b/doc/changes.rst index 1fc30c15c..b0401d6a1 100644 --- a/doc/changes.rst +++ b/doc/changes.rst @@ -5,6 +5,10 @@ desitarget Change Log 0.31.2 (unreleased) ------------------- +* Add URAT catalog information [`PR #526`_]. Includes: + * New module to retrieve URAT data from Vizier and reformat it. + * Code to match RAs/Decs to URAT, as part of that new URAT module. + * Substitute URAT PMs for GFAs where Gaia has not yet measured PMs. * Use ``MASKBITS`` instead of ``BRIGHTSTARINBLOB`` [`PR #521`_]. Also: * Extra options and checks when making and vetting bundling scripts. * Option to turn off commissioning QSO cuts to speed unit tests. @@ -15,6 +19,7 @@ desitarget Change Log .. _`PR #518`: https://github.com/desihub/desitarget/pull/518 .. _`PR #519`: https://github.com/desihub/desitarget/pull/519 .. _`PR #521`: https://github.com/desihub/desitarget/pull/521 +.. _`PR #526`: https://github.com/desihub/desitarget/pull/526 0.31.1 (2019-07-05) ------------------- From 8fe846169224328d7d4135475b025362e72a97a2 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 13:14:22 -0700 Subject: [PATCH 20/30] clean up one logging message --- bin/select_gfas | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/bin/select_gfas b/bin/select_gfas index c685a429f..dd65d60dd 100755 --- a/bin/select_gfas +++ b/bin/select_gfas @@ -80,8 +80,7 @@ gfas = select_gfas(infiles, maglim=ns.maglim, numproc=ns.numproc, tilesfile=ns.t extra = {k: v for k, v in zip(["maglim", "mindec", "mingalb"], [ns.maglim, ns.mindec, ns.mingalb])} -log.info('Writing GFAs to file...t = {:.1f} mins'.format(len(gfas), ns.dest, (time()-t0)/60.)) +log.info('Writing GFAs to file...t = {:.1f} mins'.format((time()-t0)/60.)) io.write_gfas(ns.dest, gfas, indir=ns.surveydir, indir2=ns.surveydir2, nside=nside, survey=survey, extra=extra) - log.info('{} GFAs written to {}...t = {:.1f} mins'.format(len(gfas), ns.dest, (time()-t0)/60.)) From 03b559244bf056965c22d647dd6ee49cfcc6302b Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 17:12:52 -0600 Subject: [PATCH 21/30] fix sphinx errors --- py/desitarget/io.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/py/desitarget/io.py b/py/desitarget/io.py index 770254143..baa2262ac 100644 --- a/py/desitarget/io.py +++ b/py/desitarget/io.py @@ -639,7 +639,7 @@ def write_gfas(filename, data, indir=None, indir2=None, nside=None, HEALPixels at resolution `nside`. survey : :class:`str`, optional, defaults to "?" Written to output file header as the keyword `SURVEY`. - extra :class:`dict`, optional + extra : :class:`dict`, optional If passed (and not None), write these extra dictionary keys and values to the output header. """ From c10122da8f00d38703e2ff04ae4584720f6c0a2e Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 18:08:24 -0600 Subject: [PATCH 22/30] more details about accessing the URAT data in the FORTRAN README file --- py/desitarget/fortran/README | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/py/desitarget/fortran/README b/py/desitarget/fortran/README index 9d0bafd2b..6de3f8059 100644 --- a/py/desitarget/fortran/README +++ b/py/desitarget/fortran/README @@ -12,6 +12,15 @@ v1dump: An executable compiled at NERSC via: gfortran -o v1dump ADM-v1dump.f v1sub.f -More about the UCAC/URAT survey is available at: +More about the URAT survey is available at: http://cdsarc.u-strasbg.fr/viz-bin/Cat?I/329 + + +This code is used as part of desitarget.uratmatch.make_urat_files(), +specifically, desitarget.uratmatch.urat_binary_to_csv() to download +and format URAT data from http://cdsarc.u-strasbg.fr/ftp/I/329/ + + +At NERSC, desitarget.uratmatch.make_urat_files() has already been run +and the output is in /project/projectdirs/desi/target/urat_dr1 From 40894edcecb5df5165a32d4eab0ba1d378d6c24c Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 17:37:00 -0700 Subject: [PATCH 23/30] remove fortran binary --- py/desitarget/fortran/v1dump | Bin 31064 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100755 py/desitarget/fortran/v1dump diff --git a/py/desitarget/fortran/v1dump b/py/desitarget/fortran/v1dump deleted file mode 100755 index 389e8e60a564de90cf64716480ee21689e31a92f..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 31064 zcmeHw3v^snz4w_+(x$Y{q%CLxfg=xXF`1-Kn$p4~?UV^75K03@TBg%vl8h#kkj#{p zOQ8_P+i@6(0@~%rqZiS|%JTps z>|vy=O2s(*T@~L)WEC5oaR) z6pqty6yf+Zj;T0i;9yzl&Ph1fE~w+91nrkh8%}^Q>y~vkMJbxB;^`_r39&BA$w;4q z<3t>%;$UDI^&t=+@pKdTbRKq*+fl%5_$zbpe>MmH4=^g5p0XVJ@5sSFB?rGh2mYZP z_B@qC|Hd5n({kW%$dPYT4*cRAdgvPyM`l=`gTF0@{s(gKhoCr{|KHDnkL1u_n?wJx z9QeIC^nWu4{=yviS99PE^!w(f_RdIWAk5xA-|+jJ+uFnau1KIW;`bYVP4$zdt0~y& z?+gZ7{nvGdB0(Fqp*`FXh(N3}6mGWRyCURoOvAKyQ5*4BNyT+S@t; zox$-MEc-&?NU#}pW||0H0+pGECD_ms=nOPOvWSjnu~`JY4Xy268D@vV*G5B~X~-Z8 z6A874yV9iTF@GrB7~GhK9%sL0gestoyBVs3;SHhA_HbJ;9Pw`mbcO=!TZ8@(D%7gM zE^q}po3rGOT5Gd0S~&$eI|G~8M%vgo782CRfvBOa!)Ob%wzfAAYHQyB&-lrxT3TBm6gHYVgF)k>n(C@${<-eC?ghq0 zD=uBRYGvgm{`1g&^Dyh=iNAbt))4Vb+c{D+N^nef80R3tvERfIJtc%`C*Qb}LbaYW6PyXh=ego4`o^ijDAH8%P(PI9w&to&|B%X5?<#jM^c%?#+X4&xe`K{cB z*J}g{EVAMC+)KR2hUc89gU^Q7^C)R`HoRW*5Wm`n=RB>$Ivbwzt`03Wyq#*VV z+)VsN8=mU|9k$r;TsP>j)rQw=4AOdRc)b=Pewz(%AHTQT@H5D$!b3JZ*9AH}X2Ww0 zp~Eg4p6eVPdTn^Fd35Nv;kh=_VZesxxwl zrae3$;0Flj5$+Z6w+XW;4?ia0ZxCh^9^NkCuMuX`9qtkEEri)*hqnm$dctg~!yN+d zB+Mo{yiUMhCd{TeTqof5gxMsAJp%p`VK&9#asjU-%qBQ&3ix8eY~(-~eT;5P}g$qd&C_*KGeD#IQDzeJc#WVl?w zhbm2@sRytKbN*BvfBqU@{8eA~Ya?~5s{4A43q6LfZ*R{9jQSp*=`n_W4Z(MtTxTpd zjQU-Jm#s5=OT{_<$H=F5Q4O2U)U{6&Id)j^t=%vA_Yha_yZO&bz_<4Ur!PL@+dFij z&vDTA{0EVfEDe)W8VoJJW$E8xzY+a(y)XX9p&L;F2Mdk><5;_2*t3k{&+l#&B7H+^;hBWMcJ~(gw(Vv=fHQoq<@Nl~2w>hLyW$4s^!%)x=*PT4E`T_1@_6C1N6DsJipGR*KuY0uNx^pav z+53n~CIzFTe~+mOdMgBL;%^L15oRzGWuH(Cb?*p!w7VV?Q0x&#Hg~ycqQG|d7L0e- z51?GSi}I2&k))F3rh>{R%Z0l*R9NnsYAICah>T%Qqop8OR{(eQ(voECr-+oh9$-oC zCBiCv%3X?m9e9+xE@v{Qy9A@nUG*ZLH->)mKFh&Olzl=m>^lHLpLb{l8cd(}wOgHr zanSqAOAdMumAH0#f4StUNTG-~MJM%nUtp^`G7$;HY|P9}HQ)vCiEMRfz-*O6Cm^%# zp}e^Fa9?Z<7&-Lyd0*j((-&I;e<5KcyhDk#uf=0W!JH81iv6y8dR;qvW4*4}HS9Be zvHrh~jV*-!zSsa#OCg55If(e;F_ZSkX5lQlLSM{4lA8Ns2f)k#4TV516B>$<;-K3X z>jl}U6S|mWQ=AB*=y=_m4;!v(@8S49G(q9hq~KR|@HF^f0>yUcLD`;5p_>Z6CkRQF zM$_{_S7QA(F{1Zero`G)1Zo+OS~x7D<4a{b=tAgZf2e$tWv$#LN*wQ#Lq%maU-2p| z<*O0c+ZQ{EGbQ_C4WfW*Bgs3iI;td+9P|JP{BN1UbhT&+iWzsJMAt{=is&<=6KTL_yU(;byUu+Q?IJ-+iq?WF}z-74Z*X_CljmDRJit{{5 zO4s_5^)Fy_pcH(>x~_b3D{`-nk3y%c;RmspNQjDH-xz{+*<(Tjo+_8h$T1L)m9X7p z`h?qz(ClFb_FnTl>8qT*rpIz1cOJ6#n&*$9z}};>wCoezhn?Mr9k=#ImvwJGTB?g3 zU6Ak|^=sUll_Pwv>`I1vGyUA+Kr!8_)xc5D%fzBh}hl50c#`A0fj}oQdzYP#$ z)uAP3q)M5SmVJrLzI(`#$d8#=4aB{JuATLV65c`JeA~9(RY32ETo3k)bKn)yQEcab zwS0o%Qu_q$JX`O1;y-C*QutjJLAYENa|2vwE5fe<$_CTj>xhe%GyD_SG~Sx4QW+za z5G`h0&+TX4nGNAf>5^GG{MR9-tQb;OB<|*bEs<T`hg5I)xZ*ZB68QNfoAyZyaj`)fy)@i!yO26`o~c_%<}vo*rK^ zIC6F6nzj46RKv=wIKKI~FCHE7CF(~KUL;>c+j{Puv9UcS%zBy^07Eyy;-Sx?vP8#x zMvBHaqpLzLUO(hZct_uM72|@&JA&9~a=G_zDAtgqPBRo2;~{dB>$==!nI6c;jom(0h( z?D}#22jp(D*emPhf^{ngsrXY{S=I3!TSRpCMiC`D5G{|#T1Av>L6i`HZyT0Pl~2J( zXp=+}v9IzcL|_d`>4dS=i=pDi6H_o(9j@CSAWvUx8+bhMbYE;cQfEu#AXN$EVWi~j z(t{MYLHc5MDV#4EBa(x`(zg*yaQ9izv8N8wB;=zeu?!XEP;^_H`W&e#Ef=xcA#FU( z)6gy|P3)jYWO3Enrwh8p+6Vei1P_Ljxc8NDWtj~fC^UVCA$UY+p^xIRQKGQALBxz? zi7+b=$)ZQBqRZ~YQmUm3lPp)wTh`L$6)s&~iSI)(X{uEoF@NUU9K#WN&yC74N1!0L zW2UCWn46nmW?g;Z!^@Y9-M|vXW4mBfpQx!V9$i?7a(dFGxe>-B1V5&_TVo_<=H@2S z@@6k|$`%)o{gj24j=M=X?gb!Zn_<>i^1;I!F-anU6!Y9tU^t+RjU3%eG$@-jCaf8r ziAFlkh&ePuW)-vQi0w(L&DlKf`wL7|&5xJZd*Z4Vf+H4zem+Prv*dGAnci<@3hVw+ zrn8Z0;%?C^6QT!OMHDUL&6G{d`it4Jp7&2?{pg>qV(k{?eL$5rnx5&41T)J!|DVjX z!j>ubP4?_knHuSt2C`+k;GF6I@Y(AYg%CoZU&Qz9n1xm;1qXC~tEux`Vt zn$+D}R}+7Ak11^F8}cw`DJ3so`jlh-Rp$HIG~zRg_+>$K#oovKZW|Q8LuEG~;B0|u zv~sPQtZqK~7%-SWHnIfq9UUT)tuLNj2dE|nw9chUeqlAGKV%&?_*E=HjL%50I z+h?CeBwzPviR+GAS)x8MqR2Ue3eSytbI_I$~M!z zazv0|mn^W2v{#bu(IjM{oACD}X?qq@LXvi965XI0eupGIk%d$*Nx1frX4vvtEJ?jt zNM}gWeOXBPlC)EkY-S8gZ+s^U={ZT-tw}aBek4h|vXHh((jHB+<+WLop2|X6FG&w; zk}a=`B~0M{f+ zzuYC3wb3{9i#xdyf*D<&Wh}#0Aj_7nX z;`NF9ki&j7JYcvi5o6%33*io5()%EV#cDn2#hBMOG%A*$-g}rnN37Gm+nGLDre9!s zs!aDX?U3pFnEp>GyPN4@ncl_pZ)Eymrk|JTN16UVGW`hCPs(&Z(-E0I$n?L-^dYA2 zmgxgb$7Onu>CH0z64Tep^ovXfk?!+~XVMatq2gf0wV}Se%z;7qL zo+6zBQa4l`+1Zn*e-NJJ6caB_)Nj-2x_sAai5@(6P_1!w1BEJEA60-*0)V4)O=>@~eP2WJjz9;(#5o28cmB zVjU1K*%6IEyl6+Xz=i?44Xr>tZ%1?h@e4a50>tBX#6}=~WJlZp#1CzVo-G)Qkn=0H zoO^DSxZm1wT%myPupM`s#Jy_A^+?=ncHCVO_qrXoP2z^^xO*h-O*?M8#2vNc9+bF0 z+i?#Ox6}J$d)Xe5xTo#7$AIgR8#I+y=$#D{Lnv0JgB)4oI|f8_cRwTC=0c}y!{rJSgQ6`ijlfH{@d#JUQ{}lC(J(YAV>P;6~RHP z!(XiwE>FzDy%9-nd&MGS%k&kD2d*FNuTGvg0XvldFd`&RxHVUXp@_%)=|%9uzFKYuyQi6BOB zH}Bx?8S!}F!M>rLzqbekKm98ADsMhO2aAh=%4@8BHII+`SMlg!Y>;7WtR1L{d;4qR zwY~7N*#5$nS3wi*Ka0(>3SbX&Q;gU;jXew4nSCre%I@uzGYC?}5;x3;>;Mk!F0S+# zduHN9tF93XbnoGGtWRR~PMgR~WPL6*aR2RX*EEdA(hi{(2*zaP{+<9V(|?2;TW`Cr z+B@|4?=b#i^Tj(5FIL>56yMMcU)K!%f?*(GB)1y>ap{cRfBDn=wekM<-rW2C^mzZ? z(KBZ6iw`90urcB76`cL#92@*l&SJq?oM6_238P})jZ;L5u%z|cXt+Wompv2Z-ym(F3#>oI2^MmuEaisN4dJslTUrHpchG` zm>k5`s;aa4QIQx-5eHOqI4PN75tnNJ&?1qAu^sC{r2199AHmGJ{4|Tu1sX4`CCX+} z0%q3bs{s%XA*_=lpBU@nJ9{`;Omv33h1lxsVj{OLa9oE*I8K-AtDHG1-A&%tAM4z0Sn)$@CoA1?Fas(D<$e@tO(|`3rKOAx>Qis z?0v^|GJjPC78IS5Io_abA*3~q6(EGbDr0A(OdTJK-cVMiLz>TfI0CrwLws1&9)=0 zYxYi&M$Mk9YjyzoL=E;cfttlGRX=MGmJ0l{<2utW)>1<3tKpK7s|4K9#+i?&fw>q- zD(H)`LoukH_Z6$>@kXBNd3BIU^|$W~vLZ3yhb%(?njD6K*G`yy}J%Cc+6=1&E9AgPw`yJ*f#{ zmyOJQH^5O#5&uZ|y5_X0tD?+_k15s*Xw7+3K$(L3luOn+abS0i20AFstjZb=p(F-KisS42l$_wRM> zd?xlxbaA3~)YtF~41<4{oR9!_!!tF>e3!x}H*qYXH>wg-ZjRsf)Ciu$5Leev_6T#1 z!isn;?!zYOtHdNr*?cHFPRQ|`(vX<{!!)+;VMpS&Bo$UX=eqq@Yz)fM3T^3T%TibD zM_{m9Ys(;ohUuyh)ec;+3d?1_q73o*H|v63JC-R$Yax|*Vh3$f!c{2|;Ung%d{mTD zrJ03~ggAXPyaWS54`^9jGxE`0lN=FUQdZ(ec8$s)vyim?2AX# zI9XfH&_`}LH|>UokKXmE_`|xsk?Q({ zUy>bAWs*Y6ZqEXI{I*SHvTd_~D3u;5L>G>*3lE1ezb9&k;?I3S47*>bJzDYHjpd2j z*I3#lhmzy*+ja~Ht`u$vIIhHDFh=qtZae?cEc_!D{}4v|ZcGLkT)$9zy!)>wxb8Rs zcNvk*`L)kRr-*pOb=xsKu!D5>UPrvPKN0QcSR8$~qV`>n>;8it*Un?HW3JdQFnT1S z@5ZB;OLmCXa6xK#O)*B#Ped^l9K%#F0(Z`Z2%y^JlyZr$Nv?o%cF3t<{_R!qCt}uI zfSuOp@$Nr6;?bjE5u@Egr3$kORKdr`e2({&d}7I-G&V@8LAT;P*X^GYqn)sHd18Vx zb$Nw%)Mjf<%GR_DQs|oG6z1D~(1EP09dD&&Eo*qm?P4N(5PglU`M+rQ9X_u5q=$3b zy+j$9ym2hIsfnM8SyowG6<2mXfnBV39mA)ik-_5E9Nn8=Gj1$ZZA#9iY=CHKX*f1O zw5H)7BE~WY5C7ka?uMsjudIxZRbEjYKUP^^9sgD3s+#0Erpv@n0+OYSFRz*X$2l+i z_Ws?;r}Kx}evZ<`{uDV~%h$x;tcf364jp4BzT)fN@9~OSgy_>(_>?*syJsd27VUa*9En5z|#MUEX zE&_QloFOA)!l^Wl z$A$>}HM=1YHpA@^xGLOeMp}ZxXNzHFM{6+96*TcHMM#H(*O?vo-QeO`!HuDad0i;d zVxGD2%-I{w)uym0*xA|MX)eL<8$&G51&|0fnvJ2h2J?cY=KKZXXOpJ*nc`w|{iX== zX$^Il;#Z`qNQSZH3Xe!{xuRF%7D*VutU0U8%I2&tsVrGB2UT7+du640fjM`++0q{E z>@u0RS#B;4tnV^sEnha<1hJyrYz#J=i)W#}s_`?$5{T<~#d+@evu9P#u@Dwvp?e-5 zL=!>T;koXMm4u&tIxgZZ2j_{{xaLid;4?0fvH6ICQFk)jBOB800NkM||Jm5sVZ_F>V`GI5{I(JCV*uYk zYy$qqi`bb${By)Dh~2o@;dW0m;)f7_6LCM{mwz)h#>Hb8_g;%}k^d%c1usJU6rKrK zjkp~<aW?4W_dPQ{)x2d>CtQ4+vpe5;&m`DIo-H_5{Az6Me8d*d1`r7&(1r(b ztOZYi2(2e?jbI}E0FL`Wudvf2QU~d;;rKD=Wp?_PCH-vRxKnx%bT4lkX#F)w=WY4j z1Mep(ahVO-Xb%Z?=6~jTY;~r{HqguMT@#8))skMixyQD zl~)$cswy(8ie`eps;H3Lk;W3Vy+P#E_3)kTLK+5`%u98XrGl(Cgb% zBxnrOt`jM)-xxGKUxLQQJda5I=Lci$On96U#?1v0$|p#eOG1Wqiq54f!^6s6{)(Al zri=_cnvxdOWdvSivQFTYn6X_Y^z{ONN5-Jv$pp2Kkn$BQ^F zTDEMlIqUNJ%e`g|U(Pq%obO)5mm>*0-pDuGJkLGPU2e`TpR=fZ&fEp&tYvs@lUW^Z zaL+a>YgR4yL+&g7>Tpv#)QNCKu(Jy0Hm(np2NxSx1p}Rnk#EIftvNkE zzZb8IvV{_)WtWsaOuv_s{?8UTuH+eRoOO>(=UO7@k63g)p0SP@^mvvp!EtA454|n* z^R*`oo#H(HB=h?miXi#Dsz2!QtWb@g8qUNYmy2lC2$cY-O8sY`Grie6d$i z=-5ci^AI9=4URlVVHd!>yoo_V`4d+&Dx7$^piPJhTEUkPIr9tf$}rFhZYHvD;@22C zCh&qkzY!%sG0#|uRH37r2!uPqow9Eg6&S4GLdP@|%kgtX&UE7& zU{TpS9|2F^3*3=X#^y`&3iwVz$1*L-LNPe*`2)@q_-KS<*$ROI_*dq(jFcKuej(Ps1K8_uG|lh2=AJ>$ZY^Je0Wou4zRX9#I%te=owA99+Q zsq<5oB4|xzP~di4c(Tmm^CwfzGxJ+jCp(bEgaS(<5B5z=S%Pw+CO%E3Nx&9L`^a=c z{1ykwWf+scb_@3MXMVj3CQkv2Orfp#Eu}MOmkx zhUo#*`3u95mk->xMIK+u zmPc8%t5dx=wz;7}6}7Z0+7*;1zn}e-Qr77FHB8`o~CU=W=i0<~z z{f5Iiy>PlS9C1!J<#}nTGkmerb@3XEa(U12)c1|TZ{_8|y0stI1L2eh^w{nEyyJ_) z>nOe+znx!DxN&|g(=}Pv(%u#<3#^ZZS{utcJKAKTs=m5rd08{xY*>n$19NeMDayCz znw6KEbKUdi@O_Ex;bt=$?g}-BaW8`-b5{%Az?qJs$TB(6o26}$XgIjEIT#Lhh8h3{ zIvZM+%2)K7rOjsPrE|?vICN=qI9gF*mNp??Q~_t>UA{qpfv!#Ah8DaVvOU^W8Vz5E znHNh!_ASFrR880uYSa?z(v&uV zx8W8AgtazujKVTPf zVom9BX46w(?GL2$PcZa8LOOn;k?s21iM6}FZcpc*Wa#Vvbo>d1-X}=MPqy|?((zM_ zI_1xF{8YoVt@Aa_i8rI`^?=4Pb{a10`XwD-Z0MDGI)0kC$2d3zD_#+B)m6&ogF38PQ2k=zyH=wbaI!aE*n2u^JUraCmSzh#dF;|J~*wn zwx`pa;%;O*j*wGa@n^%IJ{~Ff&&Yu{bKrGvOeZ^yGt&`(25{sXZ&~AlW^m-dPXRs; z?POg>`zccRo(w!KA^$)I{yg9v*!|FcW_}F3JaJ^H_cq|E$58u`{`@ug9eAVr ze6I(o^pH*e9|g;das7hVm=5eHW%dhMvE%#Kbm8#v{gfAqjx7DC7KcDxYJtvfSMCRgI0Hk<1J6&MXPw(ZxPtN1Jv}aoF zP)qS0RnxPKd(*+se0QGg5v<~w2c%9{`0WcMUOYPhJom-3*~9(tZ1|AUbNjhck9ZaW zdaldCpU8pd+gGybe>Dfb0RGR$@6Lg*&4CXHe39|DGo=0ESr6F%l^pzE1D^iRtoLsz zJ$IBzJ>ppn=(!(w6LDsE6!>Y_Mku^&T^i{p#3KhiPvtA>&xHOoRskFp_##7J2I`gB z@f`g8u1hw*y#;N>7~lB1VSDIN`t>}k*GuOKJw?W!W_kpBKnX1eo_h2+sdd%m$hR{G z{@;{-zDJosJTHWpzgc5@zC^X>3PClD@8{67F9-e)O3!afJwk(c_6Di9fM>mA`r+?6 z^qf(U-JdlAUu68NvXgra443ENKPLP0&s$J_o}s?&;GRw?wx{RjCi4)Og9kH_|H1%7t{2-}Aj#tEhIr>NGKB2pg( zN#o>)DfqakA0H2C4%$B_qCZBF#xB1vkzq=(F%SvByZS*?ztM&5l2Dp5^=XeZWm&8g zm-}f<8$*6RHtOeNqkevx0(NEOZ+$<5tt<8Eku)`_yguYp962Q9H&f6S@CA?@ES8(O zB_cjilFpdE@=iSgydH5I^&hrb;}YY_&E2`mQczcuKx2 zVtsxlO*F+}V+>Cl=Z_1%D?Mps6c+x(llg4AXyEe69a}jZ*Q^qYiMH4$a z3^$$&54tO>s<8>#Y`C#m>2AR$qTz1b6o#;jBb}16LEflL1N@+M23rFZP^pgAh~XBk z$Bm@BxgB66xDjWuY3J^27lVj9*rLYwmPRPikQ9}3f|Sr{MS=r1k4NgVC&L^>c=2z$&fnwo2!aRe4F+vL^I*{O5|ZmK-uu_` z>y-RDB|o4D6Tm}1X?d)WMOcP2S+%_0Z{p`A81!>05*qKnFM_~lM6|r#mwHKMv_uJN z`!%1ASA))HPBg6duVzk?g6ot#&7qyzeuXxgk)SWNyx!;9sN{7yb^cmj{h$i zh;J%+U4Qzzs#wVvt9GpUGxNVML*A_N2v$>7=v@LzoihC|_i0k?XIGU+Fz-_GTeM+H zh4z~Qjh#r)RoYCw@2J01naQtV9q)sz%O Date: Fri, 2 Aug 2019 17:46:22 -0700 Subject: [PATCH 24/30] mv fortran code to /etc/fortran; slightly weird for a resource_filename look-up --- {py/desitarget => etc}/fortran/ADM-v1dump.f | 0 {py/desitarget => etc}/fortran/README | 0 {py/desitarget => etc}/fortran/v1sub.f | 0 py/desitarget/uratmatch.py | 2 +- 4 files changed, 1 insertion(+), 1 deletion(-) rename {py/desitarget => etc}/fortran/ADM-v1dump.f (100%) rename {py/desitarget => etc}/fortran/README (100%) rename {py/desitarget => etc}/fortran/v1sub.f (100%) diff --git a/py/desitarget/fortran/ADM-v1dump.f b/etc/fortran/ADM-v1dump.f similarity index 100% rename from py/desitarget/fortran/ADM-v1dump.f rename to etc/fortran/ADM-v1dump.f diff --git a/py/desitarget/fortran/README b/etc/fortran/README similarity index 100% rename from py/desitarget/fortran/README rename to etc/fortran/README diff --git a/py/desitarget/fortran/v1sub.f b/etc/fortran/v1sub.f similarity index 100% rename from py/desitarget/fortran/v1sub.f rename to etc/fortran/v1sub.f diff --git a/py/desitarget/uratmatch.py b/py/desitarget/uratmatch.py index 8f383352f..dc7341f76 100644 --- a/py/desitarget/uratmatch.py +++ b/py/desitarget/uratmatch.py @@ -165,7 +165,7 @@ def urat_binary_to_csv(): log.info('Begin converting URAT files to CSV...t={:.1f}s' .format(time()-start)) - cmd = resource_filename('desitarget', 'fortran/v1dump') + cmd = resource_filename('desitarget', '../../etc/fortran/v1dump') os.system(cmd) log.info('Done...t={:.1f}s'.format(time()-start)) From c061c2700fb293808ac4b4048242d2fcb477136b Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 17:56:49 -0700 Subject: [PATCH 25/30] option to not add URAT if you want to circumvent it for some reason --- bin/select_gfas | 4 +++- py/desitarget/gfa.py | 26 ++++++++++++++++---------- 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/bin/select_gfas b/bin/select_gfas index dd65d60dd..df0a69760 100755 --- a/bin/select_gfas +++ b/bin/select_gfas @@ -42,6 +42,8 @@ ap.add_argument('-b', "--mingalb", type=float, default=10) ap.add_argument("--desifootprint", action='store_true', help="If sent, then limit to the current DESIMODEL Main Survey footprint (the default is to produce an 'all-sky' file)") +ap.add_argument("--nourat", action='store_true', + help="If sent, then DO NOT add URAT proper motions for Gaia sources that are missing measurable PMs") ns = ap.parse_args() @@ -74,7 +76,7 @@ if len(infiles) == 0: log.info('running on {} processors...t = {:.1f} mins'.format(ns.numproc, (time()-t0)/60.)) gfas = select_gfas(infiles, maglim=ns.maglim, numproc=ns.numproc, tilesfile=ns.tiles, - cmx=cmx, mindec=ns.mindec, mingalb=ns.mingalb) + cmx=cmx, mindec=ns.mindec, mingalb=ns.mingalb, addurat=not(ns.nourat)) # ADM extra header keywords for the output fits file. extra = {k: v for k, v in zip(["maglim", "mindec", "mingalb"], diff --git a/py/desitarget/gfa.py b/py/desitarget/gfa.py index dc625c430..dd95a6091 100644 --- a/py/desitarget/gfa.py +++ b/py/desitarget/gfa.py @@ -386,7 +386,7 @@ def _update_status(result): def select_gfas(infiles, maglim=18, numproc=4, tilesfile=None, - cmx=False, mindec=-30, mingalb=10): + cmx=False, mindec=-30, mingalb=10, addurat=True): """Create a set of GFA locations using Gaia and matching to sweeps. Parameters @@ -410,6 +410,11 @@ def select_gfas(infiles, maglim=18, numproc=4, tilesfile=None, Closest latitude to Galactic plane for output sources that do NOT match an object in the passed `infiles` (e.g. send 10 to limit to regions beyond -10o <= b < 10o)". + addurat : :class:`bool`, optional, defaults to ``True`` + If ``True`` then substitute proper motions from the URAT + catalog where Gaia is missing proper motions. Requires that + the :envvar:`URAT_DIR` is set and points to data downloaded and + formatted by, e.g., :func:`~desitarget.uratmatch.make_urat_files`. Returns ------- @@ -501,15 +506,16 @@ def _update_status(result): gfas = gfas[ind] # ADM for zero/NaN proper motion objects, add in URAT proper motions. - ii = ((np.isnan(gfas["PMRA"]) | (gfas["PMRA"] == 0)) & - (np.isnan(gfas["PMDEC"]) | (gfas["PMDEC"] == 0))) - log.info('Adding URAT for {} objects with no proper motion...t = {:.1f} mins' - .format(np.sum(ii), (time()-t0)/60)) - urat = add_urat_pms(gfas[ii], numproc=numproc) - log.info('Found an additional {} URAT objects...t = {:.1f} mins' - .format(np.sum(urat["URAT_ID"] != -1), (time()-t0)/60)) - for col in "PMRA", "PMDEC", "URAT_ID", "URAT_SEP": - gfas[col][ii] = urat[col] + if addurat: + ii = ((np.isnan(gfas["PMRA"]) | (gfas["PMRA"] == 0)) & + (np.isnan(gfas["PMDEC"]) | (gfas["PMDEC"] == 0))) + log.info('Adding URAT for {} objects with no PMs...t = {:.1f} mins' + .format(np.sum(ii), (time()-t0)/60)) + urat = add_urat_pms(gfas[ii], numproc=numproc) + log.info('Found an additional {} URAT objects...t = {:.1f} mins' + .format(np.sum(urat["URAT_ID"] != -1), (time()-t0)/60)) + for col in "PMRA", "PMDEC", "URAT_ID", "URAT_SEP": + gfas[col][ii] = urat[col] # ADM a final clean-up to remove columns that are NaN (from # ADM Gaia-matching) or that are exactly 0 (in the sweeps). From 542797aff219690bc58b408051f17dbd8e41ec60 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 17:58:58 -0700 Subject: [PATCH 26/30] add the URAT_DIR environment variable to desitarget.module --- etc/desitarget.module | 1 + 1 file changed, 1 insertion(+) diff --git a/etc/desitarget.module b/etc/desitarget.module index 9a5533e7a..37d973d86 100644 --- a/etc/desitarget.module +++ b/etc/desitarget.module @@ -77,3 +77,4 @@ setenv [string toupper $product] $PRODUCT_DIR # Add any non-standard Module code below this point. # setenv GAIA_DIR /project/projectdirs/desi/target/gaia_dr2 +setenv URAT_DIR /project/projectdirs/desi/target/urat_dr1 \ No newline at end of file From f0c24f0016f83c74a6792e3bff3552b84c67a834 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 17:59:44 -0700 Subject: [PATCH 27/30] add new line at the end of desitarget.module --- etc/desitarget.module | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/etc/desitarget.module b/etc/desitarget.module index 37d973d86..a1fbab765 100644 --- a/etc/desitarget.module +++ b/etc/desitarget.module @@ -77,4 +77,4 @@ setenv [string toupper $product] $PRODUCT_DIR # Add any non-standard Module code below this point. # setenv GAIA_DIR /project/projectdirs/desi/target/gaia_dr2 -setenv URAT_DIR /project/projectdirs/desi/target/urat_dr1 \ No newline at end of file +setenv URAT_DIR /project/projectdirs/desi/target/urat_dr1 From 42030925544a65ce729d5089395c1e9b040964f3 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 18:21:18 -0700 Subject: [PATCH 28/30] slight change to docstring regarding location of v1dump --- py/desitarget/uratmatch.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/py/desitarget/uratmatch.py b/py/desitarget/uratmatch.py index dc7341f76..ad9b6455c 100644 --- a/py/desitarget/uratmatch.py +++ b/py/desitarget/uratmatch.py @@ -155,8 +155,8 @@ def urat_binary_to_csv(): Notes ----- - The environment variable $URAT_DIR must be set. - - Relies on the executable fortran/v1dump, which is only - tested at NERSC. + - Relies on the executable etc/fortran/v1dump, which is only + tested at NERSC and might need compiled by the user. - Runs in about 40 minutes for 575 files. """ # ADM check that the URAT_DIR is set. From 42d74f900110e89a18da2f1bb5eb015c6fd23b26 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 2 Aug 2019 18:40:11 -0700 Subject: [PATCH 29/30] urat_binary_to_csv() doesn't take numproc as an input; check directory is empty before using it --- py/desitarget/uratmatch.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/py/desitarget/uratmatch.py b/py/desitarget/uratmatch.py index ad9b6455c..fe3055523 100644 --- a/py/desitarget/uratmatch.py +++ b/py/desitarget/uratmatch.py @@ -162,6 +162,14 @@ def urat_binary_to_csv(): # ADM check that the URAT_DIR is set. uratdir = _get_urat_dir() + # ADM a quick check that the csv directory is empty before writing. + csvdir = os.path.join(uratdir, 'csv') + if os.path.exists(csvdir): + if len(os.listdir(csvdir)) > 0: + msg = "{} should be empty to make URAT files!".format(csvdir) + log.critical(msg) + raise ValueError(msg) + log.info('Begin converting URAT files to CSV...t={:.1f}s' .format(time()-start)) @@ -468,7 +476,7 @@ def make_urat_files(numproc=5, download=False): log.info('Retrieved URAT files from Vizier...t={:.1f}s' .format(time()-t0)) - urat_binary_to_csv(numproc=numproc) + urat_binary_to_csv() log.info('Converted binary files to CSV...t={:.1f}s'.format(time()-t0)) urat_csv_to_fits(numproc=numproc) From eaf00cd94e660a91d01d8bfd0f52ab556df5fd75 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Tue, 6 Aug 2019 16:07:05 -0700 Subject: [PATCH 30/30] place URAT Fortran code in a more sensible place; better error logging --- {etc => py/desitarget/urat}/fortran/ADM-v1dump.f | 0 {etc => py/desitarget/urat}/fortran/README | 0 {etc => py/desitarget/urat}/fortran/v1sub.f | 0 py/desitarget/uratmatch.py | 16 ++++++++++++++-- 4 files changed, 14 insertions(+), 2 deletions(-) rename {etc => py/desitarget/urat}/fortran/ADM-v1dump.f (100%) rename {etc => py/desitarget/urat}/fortran/README (100%) rename {etc => py/desitarget/urat}/fortran/v1sub.f (100%) diff --git a/etc/fortran/ADM-v1dump.f b/py/desitarget/urat/fortran/ADM-v1dump.f similarity index 100% rename from etc/fortran/ADM-v1dump.f rename to py/desitarget/urat/fortran/ADM-v1dump.f diff --git a/etc/fortran/README b/py/desitarget/urat/fortran/README similarity index 100% rename from etc/fortran/README rename to py/desitarget/urat/fortran/README diff --git a/etc/fortran/v1sub.f b/py/desitarget/urat/fortran/v1sub.f similarity index 100% rename from etc/fortran/v1sub.f rename to py/desitarget/urat/fortran/v1sub.f diff --git a/py/desitarget/uratmatch.py b/py/desitarget/uratmatch.py index fe3055523..fe3c001d3 100644 --- a/py/desitarget/uratmatch.py +++ b/py/desitarget/uratmatch.py @@ -155,7 +155,7 @@ def urat_binary_to_csv(): Notes ----- - The environment variable $URAT_DIR must be set. - - Relies on the executable etc/fortran/v1dump, which is only + - Relies on the executable urat/fortran/v1dump, which is only tested at NERSC and might need compiled by the user. - Runs in about 40 minutes for 575 files. """ @@ -169,11 +169,23 @@ def urat_binary_to_csv(): msg = "{} should be empty to make URAT files!".format(csvdir) log.critical(msg) raise ValueError(msg) + # ADM make the directory, if needed. + else: + log.info('Making URAT directory for storing CSV files') + os.makedirs(csvdir) log.info('Begin converting URAT files to CSV...t={:.1f}s' .format(time()-start)) - cmd = resource_filename('desitarget', '../../etc/fortran/v1dump') + # ADM check the v1dump executable has been compiled. + readme = resource_filename('desitarget', 'urat/fortran/README') + cmd = resource_filename('desitarget', 'urat/fortran/v1dump') + if not (os.path.exists(cmd) and os.access(cmd, os.X_OK)): + msg = "{} must have been compiled (see {})".format(cmd, readme) + log.critical(msg) + raise ValueError(msg) + + # ADM execute v1dump. os.system(cmd) log.info('Done...t={:.1f}s'.format(time()-start))