{
"cells": [
{
"cell_type": "code",
"execution_count": 177,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"#Импортируем библиотеки, если не хватает чего-то, надо доустановить (pip install или conda )\n",
"import numpy as np \n",
"import os\n",
"import pandas as pd \n",
"import datetime\n",
"import drms\n",
"import urllib\n",
"from astropy.time import Time\n",
"import astropy.units as u\n",
"from astropy.coordinates import SkyCoord\n",
"from astropy.io import fits\n",
"import urllib, requests\n",
"import matplotlib.pyplot as plt\n",
"from sunpy.coordinates import frames\n",
"import warnings\n",
"warnings.filterwarnings(\"ignore\")\n",
"import sunpy.wcs\n",
"import sunpy.map\n",
"import pickle\n",
"import re\n",
"import matplotlib.patches as patches\n",
"import logging\n",
"from tg_tqdm import tg_tqdm\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'MDIdataset'"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"EMAIL = 'iknyazeva@gmail.com'\n",
"SAVE_PATH = 'MDIdataset'\n",
"tg_token = '831964163:AAH7SoaoqWzWIcHaS3yfdmMu-H46hhtUaXw'\n",
"tg_chat_id = 1147194\n",
"ik_chat_id = 94616973\n",
"sun_group_id = -321681009\n",
"\n",
"\n",
"def check_dataset_directory():\n",
" if not os.path.exists('MDIdataset/fragments'):\n",
" os.makedirs('MDIdataset/fragments')\n",
" \n",
" return 'MDIdataset'\n",
"check_dataset_directory()"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sunspots = pickle.load(urllib.request.urlopen('https://raw.githubusercontent.com/iknyazeva/FitsProcessing/master/sunspot_1996_2017.pkl'))\n",
"sunspots.tail(5)\n",
"len(sunspots[(sunspots.index.get_level_values(0) >= '2018-01-01')].groupby(level=0))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def request_mfits_by_date_MDI(moment, email=EMAIL, path_to_save='MDIdataset', verbose=False):\n",
" \"\"\"\n",
" Function for request fits from JSOC database\n",
" moment: pd.datetime object\n",
" return: filepath to the magnetogram\n",
" \"\"\"\n",
"\n",
" filename = 'mdi.fd_m_96m_lev182.' + moment.strftime('%Y%m%d_%H%M%S_TAI.data.fits')\n",
" filepath = os.path.join(path_to_save, filename)\n",
" \n",
" if os.path.exists(filepath):\n",
" pass\n",
" else:\n",
" c = drms.Client(email=email, verbose=verbose)\n",
" str_for_query = 'mdi.fd_M_96m_lev182' + moment.strftime('[%Y.%m.%d_%H:%M:%S_TAI]')\n",
" logger.info('Magnetogram:', str_for_query, ' will be downloaded .. ')\n",
" r = c.export(str_for_query, method='url', protocol='fits')\n",
" logger.debug(r)\n",
" r.wait()\n",
" print(r.request_url)\n",
" print(\"Download data and save to path\", filepath)\n",
" r.download(path_to_save)\n",
"\n",
" return filepath\n",
"\n",
"def read_fits_to_map(filepath, plot_show=False):\n",
" \"\"\"\n",
" read fits to sunpy object and plot in logariphmic scale\n",
" return \n",
" mymap: sunpy object\n",
" \n",
" \"\"\"\n",
" \n",
" mymap = sunpy.map.Map(filepath)\n",
" if plot_show:\n",
" plt.figure(figsize = (12,12))\n",
"\n",
"# data = np.sign(mymap.data)*np.log1p(np.abs(mymap.data))\n",
" plt.imshow(mymap.data)\n",
"# plt.imshow(data, cmap = 'gray' )\n",
" \n",
" return mymap\n",
"\n",
"def region_coord_list(datestr, sunspots_df, limit_deg=45):\n",
" \"\"\"\n",
" Function for working with sunspot_1996_2017.pkl dataframe,\n",
" return list of tuples: (datestr, NOAA number, location)\n",
" used in cropping\n",
"\n",
" args:\n",
" datestr: string for date in the format used in dataframe '2001-04-30'\n",
" sunspots_df: dataframe from file sunspot_1996_2017.pkl\n",
"\n",
" return: list of tuples\n",
" \"\"\"\n",
"\n",
" date_df = sunspots_df.loc[datestr]\n",
" date_df.index = date_df.index.droplevel()\n",
" rc_list = []\n",
" for index, row in date_df.iterrows():\n",
" if (abs(float(row.location[1:3]) <= limit_deg)) and (abs(float(row.location[4:])) <= limit_deg):\n",
" rc_list.append((pd.to_datetime(datestr, format='%Y-%m-%d'), index, row.location))\n",
"\n",
" return rc_list\n",
"\n",
"def return_pixel_from_map(mag_map, record, limit_deg=45):\n",
" '''\n",
" convert lon lat coordinate to coordinate in pixel in sun map and return it\n",
" '''\n",
"\n",
" pattern = re.compile(\"[NS]\\d{2}[EW]\\d{2}\")\n",
" assert bool(pattern.match(record)), 'Pattern should be in the same format as N20E18'\n",
" assert (abs(float(record[1:3]) <= limit_deg)) and (abs(float(record[4:])) <= limit_deg), 'Consider only regions between -{}, +{} degree'.format(limit_deg)\n",
" if record[0] == 'N':\n",
" lat = float(record[1:3])\n",
" else:\n",
" lat = -float(record[1:3])\n",
" if record[3] == 'W':\n",
" lon = float(record[4:])\n",
" else:\n",
" lon = -float(record[4:])\n",
"\n",
" hpc_coord = sunpy.wcs.convert_hg_hpc(lon, lat, b0_deg=mag_map.meta['crlt_obs'])\n",
" coord = SkyCoord(hpc_coord[0] * u.arcsec, hpc_coord[1] * u.arcsec, frame=mag_map.coordinate_frame)\n",
" pixel_pos = mag_map.world_to_pixel(coord) * u.pixel\n",
"# pixel_pos = pixel_pos.to_value()\n",
" return pixel_pos\n",
"\n",
"def crop_regions(mag_map, rc_list, delta=100, type_mag=None, plot_rec=False, plot_crop=False, limit_deg=45):\n",
" '''\n",
" Crop region by size delta and save it to disk,\n",
" if plot_rec, plot rectangle of regions on disk,\n",
" if plot_crop, plot only crop regions\n",
" '''\n",
"\n",
" # data = np.sign(mag_map.data)*np.log1p(np.abs(mag_map.data))\n",
" data = mag_map.data\n",
"\n",
" if type_mag == 'MDI':\n",
" delta = 100\n",
" if type_mag == 'HMI':\n",
" delta = 200\n",
"\n",
" if plot_rec:\n",
" fig, ax = plt.subplots(1, figsize=(12, 12))\n",
" ax.matshow(data)\n",
" plt.gray()\n",
" ax.set_title('{} magnetogram at '.format(type_mag) + rc_list[0][0].strftime('%Y-%m-%d %H:%M'))\n",
"\n",
" for record in rc_list:\n",
" pxs = return_pixel_from_map(mag_map, record[2], limit_deg).to_value()\n",
" rect = patches.Rectangle((pxs[0] - 1.25 * delta, pxs[1] - delta), 2.5 * delta, 2 * delta, linewidth=3, edgecolor='r', facecolor='none')\n",
" ax.add_patch(rect)\n",
" ax.annotate('{}.AR'.format(type_mag) + str(record[1]), xy=(pxs[0], pxs[1]), xytext=(pxs[0], pxs[1] - 50), color='yellow', fontsize='xx-large')\n",
"\n",
" plt.show()\n",
"\n",
" submaps = []\n",
" for record in rc_list:\n",
" \n",
" filename = '{}.AR{}.{}.fits'.format(type_mag, record[1], record[0].strftime('%Y-%m-%d_%H%M%S'))\n",
" print(filename)\n",
" filepath = os.path.join('{}dataset'.format(type_mag),'fragments', filename)\n",
" print(filepath)\n",
" print(os.path.exists(os.path.join('MDIdataset','fragments')))\n",
" \n",
" pxs = return_pixel_from_map(mag_map, record[2], limit_deg)\n",
" bot_l = [pxs[0] - 100*1.25*u.pixel, pxs[1] - 100*u.pixel]\n",
" top_r = [pxs[0] + 100*1.25*u.pixel, pxs[1] + 100*u.pixel]\n",
"\n",
" submap = mag_map.submap(bot_l * u.pixel, top_r * u.pixel)\n",
"\n",
" if plot_crop:\n",
" submap.peek()\n",
"\n",
" try:\n",
" submap.save(filepath)\n",
" except:\n",
" print('OUCH')\n",
" \n",
" submaps.append(submap)\n",
"\n",
" return submaps\n",
"\n",
"\n",
"def request_batch_mfits_by_date(moment,\n",
" period_of_days=30, email=EMAIL,\n",
" path_to_save='dataset',\n",
" verbose=False,\n",
" type_mag=None,\n",
" token=tg_token,\n",
" chat_id=tg_chat_id):\n",
" '''Request batch fits for a period of days and return:\n",
" request url\n",
" period of days that was apply\n",
" first date of butch\n",
" last date of batch\n",
" '''\n",
"\n",
" c = drms.Client(email=email, verbose=verbose)\n",
"\n",
" def set_str_for_query(period_of_days=period_of_days):\n",
" if type_mag == 'MDI':\n",
" str_for_query = 'mdi.fd_M_96m_lev182' + moment.strftime('[%Y.%m.%d_%H:%M:%S_TAI/{}d@24h]'.format(period_of_days))\n",
" filename_to_check = 'mdi.fd_m_96m_lev182.' + moment.strftime('%Y%m%d_%H%M%S_TAI.data.fits')\n",
" path_to_save = 'MDIdataset'\n",
" if type_mag == 'HMI':\n",
" str_for_query = 'hmi.m_720s.' + moment.strftime('[%Y.%m.%d_%H:%M:%S_TAI/{}d@24h]'.format(period_of_days))\n",
" path_to_save = 'HMIdataset'\n",
" filename_to_check = 'hmi.m_720s.' + moment.strftime('%Y%m%d_%H%M%S_TAI.magnetogram.fits')\n",
"\n",
" return str_for_query, path_to_save, filename_to_check\n",
"\n",
" str_for_query, path_to_save, filename_to_check = set_str_for_query()\n",
" print(str_for_query, path_to_save, filename_to_check)\n",
" if os.path.exists(os.path.join(path_to_save, filename_to_check)):\n",
" logger.warning('Files already exists. Skip downloads.')\n",
" return None, period_of_days, moment, moment + datetime.timedelta(days=period_of_days), period_of_days\n",
"\n",
" print('Magnetogram: {} will be downloaded ... '.format(str_for_query))\n",
"\n",
" r = c.export(str_for_query, protocol='fits')\n",
" print(r)\n",
" multiplier = 1\n",
" print(r.has_failed())\n",
" \n",
" while r.has_failed():\n",
" period_of_days -= int(10 * multiplier)\n",
" str_for_query, _, _ = set_str_for_query(period_of_days=period_of_days)\n",
" print(str_for_query)\n",
" if period_of_days < 10:\n",
" break\n",
" multiplier *= 1.5\n",
" print('Export request has failed. Reduce number of days in it')\n",
" time.sleep(2)\n",
" r = c.export(str_for_query, protocol='fits')\n",
"\n",
" print(r)\n",
" print(len(r.data))\n",
"\n",
" try:\n",
" r.wait(retries_notfound=10)\n",
" except Exception as e:\n",
" print('Can not wait anymore, skip this. Get Exception: {}'.format(e))\n",
"\n",
" print(\"Download data and save to path {}\".format(path_to_save))\n",
"\n",
" with tg_tqdm(r.urls.index, token=token, chat_id=chat_id, desc='DOWNLOAD BATCH') as batch_d:\n",
" for ind in batch_d:\n",
" try:\n",
" urllib.request.urlretrieve(r.urls.loc[ind][2], os.path.join(path_to_save, r.urls.loc[ind][1]))\n",
" except Exception as e:\n",
" print('Get error while trying download {}: {}'.format(r.urls.loc[ind][1], repr(e)))\n",
" print('Skip this file')\n",
"\n",
" first_date_batch = r.urls[0:]['record'].values[0].replace('[', ' ').split()[1].split('_')[0].replace('.', '-')\n",
" last_date_batch = r.urls[-1:]['record'].values[0].replace('[', ' ').split()[1].split('_')[0].replace('.', '-')\n",
" len_batch = len(r.urls)\n",
"\n",
" return r.request_url, period_of_days, first_date_batch, last_date_batch, len_batch"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"request_batch_mfits_by_date(new_start_moment, email=EMAIL, type_mag='MDI', path_to_save=SAVE_PATH, period_of_days=40, verbose=True)"
]
},
{
"cell_type": "code",
"execution_count": 230,
"metadata": {},
"outputs": [],
"source": [
"t = drms.Client(email='metya.tm@gmail.com', verbose=True)"
]
},
{
"cell_type": "code",
"execution_count": 235,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"( T_REC CAMERA\n",
" 0 2017.07.01_00:00:00_TAI 3\n",
" 1 2017.07.02_00:00:00_TAI 3\n",
" 2 2017.07.03_00:00:00_TAI 3\n",
" 3 2017.07.04_00:00:00_TAI 3\n",
" 4 2017.07.05_00:00:00_TAI 3\n",
" 5 2017.07.06_00:00:00_TAI 3\n",
" 6 2017.07.07_00:00:00_TAI 3\n",
" 7 2017.07.08_00:00:00_TAI 3\n",
" 8 2017.07.09_00:00:00_TAI 3\n",
" 9 2017.07.10_00:00:00_TAI 3\n",
" 10 2017.07.11_00:00:00_TAI 3\n",
" 11 2017.07.12_00:00:00_TAI 3\n",
" 12 2017.07.13_00:00:00_TAI 3\n",
" 13 2017.07.14_00:00:00_TAI 3\n",
" 14 2017.07.15_00:00:00_TAI 3\n",
" 15 2017.07.16_00:00:00_TAI 3\n",
" 16 2017.07.17_00:00:00_TAI 3\n",
" 17 2017.07.18_00:00:00_TAI 3\n",
" 18 2017.07.19_00:00:00_TAI 3\n",
" 19 2017.07.20_00:00:00_TAI 3\n",
" 20 2017.07.21_00:00:00_TAI 3\n",
" 21 2017.07.22_00:00:00_TAI 3\n",
" 22 2017.07.23_00:00:00_TAI 1\n",
" 23 2017.07.24_00:00:00_TAI 3\n",
" 24 2017.07.25_00:00:00_TAI 3\n",
" 25 2017.07.26_00:00:00_TAI 3\n",
" 26 2017.07.27_00:00:00_TAI 3\n",
" 27 2017.07.28_00:00:00_TAI 3\n",
" 28 2017.07.29_00:00:00_TAI 3\n",
" 29 2017.07.30_00:00:00_TAI 3\n",
" 30 2017.07.31_00:00:00_TAI 3\n",
" 31 2017.08.01_00:00:00_TAI 3\n",
" 32 2017.08.02_00:00:00_TAI 3\n",
" 33 2017.08.03_00:00:00_TAI 3\n",
" 34 2017.08.04_00:00:00_TAI 3\n",
" 35 2017.08.05_00:00:00_TAI 3\n",
" 36 2017.08.06_00:00:00_TAI 3\n",
" 37 2017.08.07_00:00:00_TAI 3\n",
" 38 2017.08.08_00:00:00_TAI 3\n",
" 39 2017.08.09_00:00:00_TAI 3,\n",
" magnetogram\n",
" 0 /SUM95/D993328765/S00000/magnetogram.fits\n",
" 1 /SUM94/D993330972/S00000/magnetogram.fits\n",
" 2 /SUM95/D993334611/S00000/magnetogram.fits\n",
" 3 /SUM92/D993340526/S00000/magnetogram.fits\n",
" 4 /SUM86/D993347311/S00000/magnetogram.fits\n",
" 5 /SUM89/D993353060/S00000/magnetogram.fits\n",
" 6 /SUM85/D993358897/S00000/magnetogram.fits\n",
" 7 /SUM87/D993364721/S00000/magnetogram.fits\n",
" 8 /SUM91/D993370522/S00000/magnetogram.fits\n",
" 9 /SUM90/D993376382/S00000/magnetogram.fits\n",
" 10 /SUM92/D993382237/S00000/magnetogram.fits\n",
" 11 /SUM99/D993388207/S00000/magnetogram.fits\n",
" 12 /SUM99/D993394237/S00000/magnetogram.fits\n",
" 13 /SUM85/D993396875/S00000/magnetogram.fits\n",
" 14 /SUM88/D993401105/S00000/magnetogram.fits\n",
" 15 /SUM93/D993403292/S00000/magnetogram.fits\n",
" 16 /SUM87/D993405560/S00000/magnetogram.fits\n",
" 17 /SUM84/D993407885/S00000/magnetogram.fits\n",
" 18 /SUM87/D993410117/S00000/magnetogram.fits\n",
" 19 /SUM96/D993414609/S00000/magnetogram.fits\n",
" 20 /SUM87/D993416886/S00000/magnetogram.fits\n",
" 21 /SUM84/D993419197/S00000/magnetogram.fits\n",
" 22 /SUM97/D993421386/S00000/\n",
" 23 /SUM88/D993424630/S00000/magnetogram.fits\n",
" 24 /SUM91/D993426960/S00000/magnetogram.fits\n",
" 25 /SUM84/D993927975/S00000/magnetogram.fits\n",
" 26 /SUM92/D993930279/S00000/magnetogram.fits\n",
" 27 /SUM85/D993932578/S00000/magnetogram.fits\n",
" 28 /SUM97/D993934956/S00000/magnetogram.fits\n",
" 29 /SUM85/D993937269/S00000/magnetogram.fits\n",
" 30 /SUM97/D993939721/S00000/magnetogram.fits\n",
" 31 /SUM97/D979843046/S00000/magnetogram.fits\n",
" 32 /SUM88/D979856508/S00000/magnetogram.fits\n",
" 33 /SUM87/D982548821/S00000/magnetogram.fits\n",
" 34 /SUM90/D982546370/S00000/magnetogram.fits\n",
" 35 /SUM92/D982548964/S00000/magnetogram.fits\n",
" 36 /SUM91/D982548765/S00000/magnetogram.fits\n",
" 37 /SUM86/D982548719/S00000/magnetogram.fits\n",
" 38 /SUM95/D982547153/S00000/magnetogram.fits\n",
" 39 /SUM84/D982549491/S00000/magnetogram.fits)"
]
},
"execution_count": 235,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"str_for_queryHMI40d = 'hmi.m_720s[2017.07.01_00:00:00_TAI/40d@24h]'\n",
"str_for_queryHMI = 'hmi.m_720s[2011.04.04_00:00:00_TAI]'\n",
"# t.info('hmi.m_720s')\n",
"t.query(str_for_queryHMI40d, pkeys=1, seg='magnetogram')"
]
},
{
"cell_type": "code",
"execution_count": 236,
"metadata": {},
"outputs": [],
"source": [
"str_for_query = 'mdi.fd_M_96m_lev182[2005.04.04_00:00:00_TAI]'\n",
"str_for_query2d = 'mdi.fd_M_96m_lev182[2005.04.04_00:00:00_TAI/2d@24h]'\n",
"str_for_query30d = 'mdi.fd_M_96m_lev182[2005.04.04_00:00:00_TAI/30d@24h]'\n",
"str_for_query100d = 'mdi.fd_M_96m_lev182[2005.04.04_00:00:00_TAI/100d@24h]'\n",
"str_for_query300d = 'mdi.fd_M_96m_lev182[2005.04.04_00:00:00_TAI/300d@24h]'\n",
"str_for_query330d = 'mdi.fd_M_96m_lev182[2005.04.04_00:00:00_TAI/330d@24h]'\n",
"str_for_query365d = 'mdi.fd_M_96m_lev182[2005.04.04_00:00:00_TAI/365d@24h]'\n",
"str_for_query400d = 'mdi.fd_M_96m_lev182[2005.04.04_00:00:00_TAI/400d@24h]'\n",
"# str_for_queryHMI30d = 'hmi.M_720s[2005.04.04_00:00:00_TAI/400d@24h]'\n",
"# trurl = t.export(str_for_query2d, method='url')\n",
"# trquick = t.export(str_for_query30d)\n",
"# trfits = t.export(str_for_query30d, protocol='fits')\n",
"# trquickfits = t.export(str_for_query2d, method='url_quick', protocol='fits')\n",
"# trurlfits = t.export(str_for_query300d, protocol='fits')\n",
"trhmi = t.export(str_for_queryHMI40d, protocol='fits')"
]
},
{
"cell_type": "code",
"execution_count": 293,
"metadata": {},
"outputs": [],
"source": [
"# import urllib\n",
"trhmi.wait(sleep=10, retries_notfound=10)\n",
"trhmi.request_url\n",
"# urllib.request.urlretrieve(trurlfits.urls.url[5], trurlfits.urls.filename[5])\n",
"trhmi.has_failed()\n",
"'.'.join(trhmi.urls.filename[2].split('.')[:3] + trhmi.urls.filename[2].split('.')[4:])\n",
"trhmi.urls.loc[2][1]\n",
"# trhmi.urls.url[4]\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"with fits.open(trurlfits.urls.filename[5]) as fit_file:\n",
" print(fit_file[1].header)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tm = sunpy.map.Map(trurlfits.urls.filename[5])\n",
"tm.fits_header"
]
},
{
"cell_type": "code",
"execution_count": 1270,
"metadata": {},
"outputs": [],
"source": [
"str_for_query1m = 'mdi.fd_M_96m_lev182[2005.04.04_00:00:00_TAI/3y]'\n",
"# t.query('mdi.fd_M_96m_lev182[1998.06.25_00:00:00_TAI/300d@24h]', pkeys=1)\n",
"# str_for_query\n",
"r = t.export('mdi.fd_M_96m_lev182[1998.01.27_00:00:00_TAI/10d@24h]')"
]
},
{
"cell_type": "code",
"execution_count": 1271,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 1271,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"r.status"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rc_list = region_coord_list('1996-06-04', sunspots, limit_deg=45)\n",
"rc_list"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"filepath = request_mfits_by_date_MDI(rc_list[0][0], email = EMAIL, path_to_save = 'MDIdataset', verbose=True)\n",
"filepath"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mymap = read_fits_to_map(filepath, plot_show=False)\n",
"crop_regions(mymap, rc_list, plot_rec=False, limit_deg=45, type_mag='MDI', plot_crop=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def pr(value):\n",
" import inspect\n",
" callers_local_vars = inspect.currentframe().f_back.f_locals.items()\n",
" print([var_name for var_name, var_val in callers_local_vars if var_val is value][0], value)\n",
"\n",
"rec = rc_list[0]\n",
"\n",
"plt.figure(figsize=(20,20))\n",
" \n",
"lon, lat = lonlat(mymap, rec[2], limit_deg=45)\n",
"print('lon, lat', lon, lat)\n",
"\n",
"hpc_coord = sunpy.wcs.convert_hg_hpc(lon, lat, b0_deg=mymap.meta['crlt_obs'])\n",
"print('hpc_coord', hpc_coord)\n",
"coord = SkyCoord(hpc_coord[0] * u.arcsec, hpc_coord[1] * u.arcsec, frame=mymap.coordinate_frame)\n",
"print('coord', coord)\n",
"pixel_pos_u = mymap.world_to_pixel(coord) * u.pixel\n",
"print('pixel_pos_u', pixel_pos_u)\n",
"pixel_pos_u_q = mymap.world_to_pixel(coord)\n",
"pr(pixel_pos_u_q)\n",
"pop1 = [pixel_pos_u_q[0] - 100*1.25*u.pixel, pixel_pos_u_q[1] - 100*u.pixel]\n",
"pop2 = [pixel_pos_u_q[0] + 100*1.25*u.pixel, pixel_pos_u_q[1] + 100*u.pixel]\n",
"pr(pop1)\n",
"pr(pop2)\n",
"pixel_pos = pixel_pos_u.to_value()\n",
"print('pixel_pos', pixel_pos)\n",
"bot_l = [pixel_pos[0] - 100 * 1.25, pixel_pos[1] - 100]\n",
"print('bot_l', bot_l)\n",
"top_r = [pixel_pos[0] + 100 * 1.25, pixel_pos[1] + 100]\n",
"print('top_r', top_r)\n",
"\n",
"\n",
"pxl_wrd = mymap.pixel_to_world(pixel_pos[0]*u.pixel, pixel_pos[1]*u.pixel)\n",
"print('\\npxl_wrd', pxl_wrd)\n",
"\n",
"bot_l_w = mymap.pixel_to_world(bot_l[0]*u.pixel, bot_l[1]*u.pixel)\n",
"pr(bot_l_w)\n",
"top_r_w = mymap.pixel_to_world(top_r[0]*u.pixel, top_r[1]*u.pixel)\n",
"pr(top_r_w)\n",
"\n",
"pr(rec)\n",
"mymap.submap(pop1*u.pixel, pop2*u.pixel).peek()\n",
"# mymap.submap(pop1*u.pixel, pop2*u.pixel).save('MDIdataset/fragments/pidor2.fits')\n",
"mymap.submap(pop1*u.pixel, pop2*u.pixel)\n",
"\n",
"# mymap.save()\n",
"# mymap.submap(bot_l*u.pixel, top_r*u.pixel).plot()\n",
"# mymap.submap(bot_l*u.pixel, top_r*u.pixel)\n",
"\n",
"# mymap.submap(bot_l_w, top_r_w).peek()\n",
"# mymap.submap(bot_l_w, top_r_w)\n",
"# hg_coord = sunpy.wcs.convert_hpc_hg(pixel_pos[0]*u, pixel_pos[1], b0_deg=mymap.meta['crlt_obs']\n",
"# )\n",
"# print(hg_coord)\n",
"# revert_coord = SkyCoord(hg_coord[0], pixel_pos[1]*u.pixel)\n",
"\n",
"\n",
"# pxl = return_pixel_from_map(mymap, rc_list[0][2])\n",
"# (pxl[0] - 1.25 * 100, pxl[1] - 100)\n",
"# # (, 2.5 * delta, 2 * delta,)\n",
"# mymap.pixel_to_world(pxl[0]*u.pixel, pxl[1]*u.pixel)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def lonlat(mag_map, record, limit_deg=45):\n",
" pattern = re.compile(\"[NS]\\d{2}[EW]\\d{2}\")\n",
" assert bool(pattern.match(record)), 'Pattern should be in the same format as N20E18'\n",
" assert (abs(float(record[1:3]) <= limit_deg)) and (abs(float(record[4:])) <= limit_deg), 'Consider only regions between -{45}, +{45} degree'.format(limit_deg)\n",
" if record[0] == 'N':\n",
" lat = float(record[1:3])\n",
" else:\n",
" lat = -float(record[1:3])\n",
" if record[3] == 'W':\n",
" lon = float(record[4:])\n",
" else:\n",
" lon = -float(record[4:])\n",
" \n",
" return lon, lat"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"datestr = '2003-04-10'\n",
"moment = pd.to_datetime(datestr, format='%Y-%m-%d')\n",
"filepath=os.path.join(SAVE_PATH, 'mdi.fd_m_96m_lev182.' + moment.strftime('%Y%m%d_%H%M%S_TAI.data.fits'))\n",
"mymap = read_fits_to_map(filepath, plot_show=False)\n",
"rc_list = region_coord_list(datestr, limit_deg=60)\n",
"crop_regions(mymap, rc_list)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rec[0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"with fits.open('MDIdataset/fragments/pidor.fits') as gon:\n",
" plt.imshow(gon[0].data)\n",
"\n",
"gon2 = sunpy.map.Map('MDIdataset/fragments/pidor2.fits')\n",
"gon2.peek()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import time\n",
"import tqdm\n",
"def date_compare(date):\n",
" return date < datetime.datetime.fromtimestamp(time.mktime(time.strptime('2010-05-01', '%Y-%m-%d')))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for date, df in tqdm.tqdm_notebook(sunspots['1996-06-04':'1996-06-08'].groupby(level=0), desc=\"TTTTT\", postfix='RRRRR'):\n",
" rc = region_coord_list(str(date), df)\n",
" if not rc:\n",
" continue\n",
" if date_compare(date):\n",
" filepath = request_mfits_by_date_MDI(rc[0][0], email = EMAIL, path_to_save = 'MDIdataset', verbose=True)\n",
" print('\\n\\n', filepath)\n",
" sun_map = read_fits_to_map(filepath, plot_show=True)\n",
" crop_regions(sun_map, rc, plot_rec=True, plot_crop=True, type_mag='MDI')\n",
" else:\n",
" filepath = request_mfits_by_date_HMI(rc[0][0], email = EMAIL, path_to_save = 'HMIdataset', verbose=True)\n",
" sun_map = read_fits_to_map(filepath, plot_show=True)\n",
" crop_regions(sun_map, rc, plot_rec=True, plot_crop=True, type_mag='HMI')\n",
" \n",
"# os.remove(filepath)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sunspots['1996-06-04':'1996-06-20']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for row in tqdm.tqdm_notebook(sunspots.loc['1996-06-04':'1996-06-08'].iterrows(), total=len(sunspots.loc['1996-06-04':'1996-06-08'])):\n",
" print(row)\n",
" time.sleep(1)"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"33751\n"
]
},
{
"data": {
"text/plain": [
"2507"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# first_date_batch = trurlfits.urls[0:]['record'].values[0].replace('[', ' ').split()[1].split('_')[0].replace('.', '-')\n",
"# last_date_batch = trurlfits.urls[-1:]['record'].values[0].replace('[', ' ').split()[1].split('_')[0].replace('.', '-')\n",
"# print(first_date_batch)\n",
"# print(last_date_batch)\n",
"\n",
"# sunspots[sunspots.index[0]>'2011']\n",
"# sunspots[date>lastrecord]\n",
"# for date, df in sunspots.loc['2011'].groupby(level=0):\n",
"# date, df\n",
"# sunspots.loc['2006-01-30':'2010-01-01']\n",
"# sunspots[(sunspots.index.get_level_values(0) >= firstrecorddate) & (sunspots.index.get_level_values(0) <= lastrecorddate)]\n",
"start_moment = sunspots.index.get_level_values(0)[0]\n",
"is_end = len(sunspots[(sunspots.index.get_level_values(0) >= start_moment)])\n",
"print(is_end)\n",
"new_start_moment = start_moment + datetime.timedelta(days=5230)\n",
"moment = sunspots[(sunspots.index.get_level_values(0) > new_start_moment)].index.get_level_values(0)[0]\n",
"moment\n",
"sunspots[(sunspots.index.get_level_values(0) >= moment)].groupby(level = 0)\n",
"# len(sunspots.groupby(level = 0))"
]
},
{
"cell_type": "code",
"execution_count": 1025,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1996-06-04 00:00:00\n"
]
}
],
"source": [
"moment = sunspots.index.get_level_values(0)[0]\n",
"period_of_days = 300\n",
"str_for_query = 'hmi.m_720s.' + moment.strftime('[%Y.%m.%d_%H:%M:%S_TAI/{}d@24h]'.format(period_of_days))\n",
"str_for_query\n",
"\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 960,
"metadata": {},
"outputs": [],
"source": [
"def determine_type_str_query(period_of_days=period_of_days):\n",
" str_for_query = 'mdi.fd_M_96m_lev182' + moment.strftime('[%Y.%m.%d_%H:%M:%S_TAI/{}d@24h]'.format(period_of_days))\n",
" return str_for_query"
]
},
{
"cell_type": "code",
"execution_count": 966,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Export request has failed. Reduce number of days in it\n"
]
}
],
"source": [
"period_of_days = 365\n",
"str_for_query = determine_type_str_query(period_of_days)\n",
"multiplier = 1\n",
"while tr.has_failed():\n",
" print('Export request has failed. Reduce number of days in it')\n",
" period_of_days -= 10 * multiplier\n",
" multiplier *= 2\n",
" str_for_query = determine_type_str_query(period_of_days=period_of_days)\n",
" tr = t.export(str_for_query, protocol='fits')"
]
},
{
"cell_type": "code",
"execution_count": 1017,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'mdi.fd_m_96m_lev182.19960614_000000_TAI.data.fits'"
]
},
"execution_count": 1017,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tr.urls.loc[10][1]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 1116,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"00:10 in total: 100%|0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000| 10/10 [00:10<00:00, 1.01s/it]\n"
]
}
],
"source": []
},
{
"cell_type": "code",
"execution_count": 410,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" | \n",
" number_of_ss | \n",
" size | \n",
" class | \n",
" location | \n",
" magn_class | \n",
"
\n",
" \n",
" | date | \n",
" region | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 2010-07-04 | \n",
" 11084 | \n",
" 1 | \n",
" 150 | \n",
" HSX | \n",
" S19W19 | \n",
" A | \n",
"
\n",
" \n",
" | 2010-07-05 | \n",
" 11084 | \n",
" 1 | \n",
" 110 | \n",
" HSX | \n",
" S19W32 | \n",
" A | \n",
"
\n",
" \n",
" | 2010-07-06 | \n",
" 11084 | \n",
" 1 | \n",
" 100 | \n",
" HSX | \n",
" S19W46 | \n",
" A | \n",
"
\n",
" \n",
" | 11086 | \n",
" 2 | \n",
" 10 | \n",
" BXO | \n",
" N18W53 | \n",
" B | \n",
"
\n",
" \n",
" | 2010-07-07 | \n",
" 11084 | \n",
" 1 | \n",
" 110 | \n",
" HSX | \n",
" S19W59 | \n",
" A | \n",
"
\n",
" \n",
"
\n",
"