Magnetic pulse of cities

↑ Click on “Open in Colab” to execute this notebook live on Google Colab.

In this notebook, we will go over the entire analysis of the data acquired from our magnetometers. Any questions regarding this project and the data used here can be redirected to the corresponding author, Vincent Dumont .

Getting Started

Data public access

The totally of the data used in this project is currently stored on a shared Google Drive folder. Click here to access the shared data repository.

After clicking the above link, the shared folder (called CityMagData ) will be listed under the Shared with me section of your Google Drive, you can then do a right click on that folder and click Add shortcut to Drive , this will add the folder to your My Drive repository.

Finally, executing the following cell will mount your Google Drive to this instance of Google Colab and give you access to the data:

from google.colab import drive
drive.mount('/content/drive')

Mounted at /content/drive

Software installation

The following commands can then be run to install the mlpy and nuri Python libraries on the Google server:

%%capture
!git clone https://gitlab.com/citymag/analysis/nuri
!sh nuri/scripts/install_on_colab.sh
import nuri

Data downsampling

To ease the reproduction of the results, the downsampled data have already been made available (see beginning of this notebook to access the 1Hz downsampled data).

Part of this analysis was performed on data resampled down to 1 Hz. We used our NURI software to do the downsampling. The code below will look through all available data within given period and save 1Hz data files to a destination folder of your choice. In order to avoid Google Colab to crash because of memory usage, we reset all the variables after each iteration.

# input_dir = '/content/drive/MyDrive/CityMagData/3960Hz/'
# nuri.bulk_decimate(input_dir,'1Hz/',station=2,tmin='2016-3-14',tmax='2016-4-11')

ATTENTION : The above command will take several hours (even probably more than a day) to complete. You can change the range of time, i.e., tmin and/or tmax , and uncomment the lines to test it.

Magnetometer data

We will start the analysis by taking a look at the overall decimated time series for both cities. The data were downsampled at 1Hz. The magnetic field data can be extracted using the get_data function available from the nuri package as follows:

Berkeley data

import time
t0 = time.time()
berkeley = nuri.get_data('2016-3-14','2016-4-11','/content/drive/MyDrive/CityMagData/1Hz',station=2,sensor='biomed',scalar=True)
print('Berkeley data loaded in %.2f seconds'%(time.time()-t0))

Berkeley data loaded in 26.70 seconds

Brooklyn data

import time
t0 = time.time()
brooklyn = nuri.get_data('2018-5-7','2018-6-4','/content/drive/MyDrive/CityMagData/1Hz',station=4,sensor='biomed',scalar=True)
print('Brooklyn data loaded in %.2f seconds'%(time.time()-t0))

Brooklyn data loaded in 25.20 seconds

Geomagnetic field

We use the Geomag Algorithms to retrieve the geomagnetic field from the USGS stations. The following pip command should usually be enough to install the program:

%%capture
!pip install git+https://github.com/usgs/geomag-algorithms.git

USGS/Fresno data

USGS/Fresno(FRN) station at Fresno, CA, 3 hours (200 miles) south east of Berkeley.

import sys
import geomagio
from obspy.core import UTCDateTime
input_factory = geomagio.edge.EdgeFactory()
timeseries = input_factory.get_timeseries(
    observatory = 'FRN',
    channels = ('H', 'E', 'Z', 'F'),
    type = 'variation',
    interval = 'minute',
    starttime = UTCDateTime('2016-03-14T00:00:00Z'),
    endtime = UTCDateTime('2016-04-10T23:59:59Z'))
print(timeseries)
frn = timeseries[3].data

USGS/Fredericksburg data

USGS/Fredericksburg(FRD) station at Corbin, VA, 5 hours (300 miles) south west of CUSP.

import sys
import geomagio
from obspy.core import UTCDateTime
input_factory = geomagio.edge.EdgeFactory()
timeseries = input_factory.get_timeseries(
    observatory = 'FRD',
    channels = ('H', 'E', 'Z', 'F'),
    type = 'variation',
    interval = 'minute',
    starttime = UTCDateTime('2018-05-07T00:00:00Z'),
    endtime = UTCDateTime('2018-06-03T23:59:59Z'))
print(timeseries)
frd = timeseries[3].data

4 Trace(s) in Stream:
NT.FRD.R0.H | 2018-05-07T00:00:00.000000Z - 2018-06-03T23:59:00.000000Z | 60.0 s, 40320 samples
NT.FRD.R0.E | 2018-05-07T00:00:00.000000Z - 2018-06-03T23:59:00.000000Z | 60.0 s, 40320 samples
NT.FRD.R0.Z | 2018-05-07T00:00:00.000000Z - 2018-06-03T23:59:00.000000Z | 60.0 s, 40320 samples
NT.FRD.R0.F | 2018-05-07T00:00:00.000000Z - 2018-06-03T23:59:00.000000Z | 60.0 s, 40320 samples

Comparing periods

Here, we plot and compare the geomagnetic data measured from both USGS stations during both periods. First, we need to extract the geomagnetic data from the Fresno/USGS station during the period when the sensors were in Brooklyn, that is, from May 7 to June 3, 2018, as well as the geomagnetic data from the Fredericksburg/USGS station during the period when the sensors were in Berkeley, that is, from March 14 to April 10, 2016.

import sys
import geomagio
from obspy.core import UTCDateTime
input_factory = geomagio.edge.EdgeFactory()
frn_period2 = input_factory.get_timeseries(
    observatory = 'FRN',
    channels = ('H', 'E', 'Z', 'F'),
    type = 'variation',
    interval = 'minute',
    starttime = UTCDateTime('2018-05-07T00:00:00Z'),
    endtime = UTCDateTime('2018-06-03T23:59:59Z'))[3].data
frd_period1 = input_factory.get_timeseries(
    observatory = 'FRD',
    channels = ('H', 'E', 'Z', 'F'),
    type = 'variation',
    interval = 'minute',
    starttime = UTCDateTime('2016-03-14T00:00:00Z'),
    endtime = UTCDateTime('2016-04-10T23:59:59Z'))[3].data

Below we show the results. We notice a global decrease of the geomagnetic field of about 0.2 μT from March 2016 to May 2018. The change in the Earth’s magnetic field is the same for both stations, located on the West and East Coast of the USA. Given that urban magnetic field variations occur at an order of magnitude higher, the variations of the natural magnetic background and its differences between the two locations are therefore negligibly small on the scale of the urban variations studied in this work.

import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from datetime import datetime, timedelta
plt.style.use('seaborn')
fig,ax = plt.subplots(3,2,figsize=(15,6),dpi=100,sharey='row',sharex='col',gridspec_kw={'height_ratios': [1, 2, 2]})
periods = [pd.date_range(start='2016-03-14',end='2016-04-11',periods=len(frn)),
           pd.date_range(start='2018-05-07',end='2018-06-04',periods=len(frd))]
for i,(ts1,ts2) in enumerate([[frn,frn_period2],[frd_period1,frd]]):
  ax[i+1][0].plot(periods[0],ts1/1000,color='tomato')
  ax[i+1][1].plot(periods[1],ts2/1000,color='tomato')
  ax[i+1][0].set_ylabel(r'Geomagnetic Field [$\mathrm{\mu T}$]'+'\n'+('FRN' if i==0 else 'FRD')+'/USGS station')
  ax[i+1][0].set_xlim(periods[0][0],periods[0][-1])
  ax[i+1][1].set_xlim(periods[1][0],periods[1][-1])
  ax[i+1][0].xaxis.set_major_formatter(mdates.DateFormatter('%A\n%Y-%m-%d'))
  ax[i+1][1].xaxis.set_major_formatter(mdates.DateFormatter('%A\n%Y-%m-%d'))
ax[0][0].plot(periods[0],(frd_period1-frn)/1000,color='navy',lw=1)
ax[0][1].plot(periods[1],(frd-frn_period2)/1000,color='navy',lw=1)
ax[0][0].set_ylabel(r'Difference [$\mathrm{\mu T}$]'+'\nFRD - FRN')
plt.tight_layout()
plt.savefig('usgs.pdf')
plt.show()

Fig. 1: Time series analysis

The following function will be run to plot both datasets in separate subplots:

import numpy
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from datetime import datetime, timedelta
fig = plt.figure(figsize=(12,6),dpi=300)
ax1 = fig.add_axes([0.10,0.68,0.69,0.30])
ax2 = fig.add_axes([0.10,0.58,0.69,0.08],sharex=ax1)
ax3 = fig.add_axes([0.80,0.58,0.15,0.40])
ax4 = fig.add_axes([0.10,0.19,0.69,0.30])
ax5 = fig.add_axes([0.10,0.09,0.69,0.08],sharex=ax4)
ax6 = fig.add_axes([0.80,0.09,0.15,0.40],sharex=ax3,sharey=ax3)
ax = [[ax1,ax2,ax3],[ax4,ax5,ax6]]
for n,(data,city,usgs) in enumerate([[berkeley,'Berkeley',frn],[brooklyn,'Brooklyn',frd]]):
  # Prepare times and rebin data
  times = [datetime.fromtimestamp(tstamp) for tstamp in data.times.value]
  data_small = numpy.empty((0,2))
  for i in range(0,len(data),3600):
    data_small = numpy.vstack((data_small,[times[i],numpy.average(data[i:i+3600])]))
  # Plot magnetic field
  ax[n][0].plot(times,data,lw=0.07)
  ax[n][0].plot(data_small[:,0],data_small[:,1],lw=1,color='yellow')
  ax[n][0].set_xlim(times[0],times[-1])
  ax[n][0].set_ylim([48.5,51.5] if n==0 else [87.5,97.5])
  # Highlight weekends
  for t0 in numpy.arange(times[0]+timedelta(days=5),times[-1],timedelta(days=7)):
    t0 = t0.astype('M8[ms]').astype('O')
    ax[n][0].axvspan(t0,t0+timedelta(days=2),alpha=0.1,color='lime')
  # Highlight Memorial day
  if city=='Brooklyn':
    ax[n][0].axvspan(times[0]+timedelta(days=21),times[0]+timedelta(days=22),alpha=0.1,color='red')
  # Plot geomagnetic field
  ax[n][1].plot(numpy.arange(times[0],times[-1],timedelta(minutes=1)),usgs/1000,lw=1,color='orange')
  ax[n][1].set_xlim(times[0],times[-1])
  if n==1:
    ax[n][1].set_ylim([50.9,51.02])
  # Setting up tickmarks and labels
  ax[n][1].xaxis.set_major_formatter(mdates.DateFormatter('%A\n%Y-%m-%d'))
  ax[n][1].xaxis.set_ticks(numpy.arange(times[0],times[-1],timedelta(days=7)))
  ax[n][2].hist(data-numpy.mean(data),bins=30,range=[-6,6],alpha=0.8,edgecolor='navy')
  ax[n][2].hist(usgs/1000-numpy.median(usgs/1000),bins=200,range=[-6,6],color='orange')
  ax[n][2].set_yscale('log', nonposy='clip')
  ax[n][2].yaxis.tick_right()
# Final fix
# ax[0][0].xaxis.tick_top()
ax[1][2].set_xlabel(r'Standard Deviation [$\mathrm{\mu T}$]',fontsize=12)
plt.setp(ax[0][0].get_xticklabels(), visible=False)
plt.setp(ax[0][2].get_xticklabels(), visible=False)
plt.setp(ax[1][0].get_xticklabels(), visible=False)
fig.text(0.02,0.78,r'Magnetic Field [$\mathrm{\mu T}$]',rotation=90,ha='center',va='center')
fig.text(0.05,0.83,'Berkeley',rotation=90,ha='center',va='center')
fig.text(0.05,0.62,'FRN',rotation=90,ha='center',va='center')
fig.text(0.02,0.29,r'Magnetic Field [$\mathrm{\mu T}$]',rotation=90,ha='center',va='center')
fig.text(0.05,0.34,'Brooklyn',rotation=90,ha='center',va='center')
fig.text(0.05,0.13,'FRD',rotation=90,ha='center',va='center')
plt.savefig('timeseries.png')
plt.show()

Fig. 2: Variances

dict = {'Berkeley':{'data':berkeley.value,
                    'xlim':[48.8,50.3],
                    'mode': 'indskew',
                    'bins': 45,
                    'daytime': [[0,1],[4.5,24]]
                    },
        'Brooklyn':{'data':brooklyn.value,
                    'xlim':[90,94.5],
                    'mode': 'indskew',
                    'bins': 45,
                    'daytime': [[7,23]]
                    }
        }

import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
fig,ax = plt.subplots(2,2,figsize=(12,6),dpi=300,sharey='row')
for i,name in enumerate(dict.keys()):
  for period,color in [['weekday','navy'],['weekend','gold']]:
    times,var = nuri.get_variance(dict[name]['data'],period)
    ax[0][i].plot(times,var,color=color,label=period)
    ax[0][i].fill_between(times,var,color=color,alpha=0.2)
  if i==0:
    ax[0][i].set_ylabel(r'Daily Variance [$\mathrm{\mu T}^2$]',labelpad=10)
    ax[0][i].axvline(1,ls='dashed',lw=2,color='black',alpha=0.8)
    ax[0][i].axvline(4.5,ls='dashed',lw=2,color='black',alpha=0.8)
  else:
    ax[0][i].axvline(23,ls='dashed',lw=2,color='black',alpha=0.8)
    ax[0][i].axvline(7,ls='dashed',lw=2,color='black',alpha=0.8)
  ax[0][i].set_yscale('log')
  ax[0][i].set_xlim([0,24])
  # ax[0][i].set_ylim(ymin=0)
  ax[0][i].set_xlabel(r'%s Time [hour]'%name,labelpad=10)
  ax[0][i].legend(loc=(0.55,0.05),handlelength=1,handleheight=0.07)
  # ax[0][i].xaxis.set_label_position('top')
  # ax[0][i].xaxis.tick_top()
  # if i==1:
  #   ax[0][i].yaxis.tick_right()
for i,name in enumerate(dict.keys()):
  bins = dict[name]['bins']
  data = dict[name]['data']
  xlim = dict[name]['xlim']
  mode = dict[name]['mode']
  daytime = dict[name]['daytime']
  day, night = nuri.day_and_night(data,daytime)
  x,y,y1,y2 = nuri.gaussian_fit(mode,data,bins,xlim,reference=[day,night],order=2)
  ax[1][i].hist(data,bins,range=xlim,color='navy',alpha=0.333,edgecolor='navy')
  ax[1][i].hist(day,bins,range=xlim,color='orangered',alpha=0.333)
  ax[1][i].hist(night,bins,range=xlim,color='turquoise',alpha=0.333)
  ax[1][i].plot(x, y1, lw=2, color='orangered',label='Daytime')
  ax[1][i].plot(x, y2, lw=2, color='turquoise',label='Nighttime')
  ax[1][i].plot(x, y, lw=2, color='navy',label='Full day')
  ax[1][i].get_yaxis().set_major_formatter(ticker.FuncFormatter(lambda x, p: '{}k'.format(int(x/1e3))))
  ax[1][i].set_xlabel(r'Magnetic Field [$\mathrm{\mu}$T]')
  ax[1][i].legend(loc='best',handlelength=1,handleheight=0.07)
  # if i==1:
  #   ax[1][i].yaxis.tick_right()
  if i==0:
    ax[1][i].set_ylabel('Number of samples',labelpad=10)
  ax[1][i].set_ylim(ymax=150000)
plt.tight_layout()
plt.savefig('vardist.png')
plt.show()

1/1 Gaussian (skewed) : amp=67980±407, mu=49.361±0.006, sigma=0.281±0.005, skewness=1.39±0.08
1/1 Gaussian (skewed) : amp=11488±654, mu=49.925±0.009, sigma=0.212±0.018, skewness=-4.61±1.48
1/1 Gaussian (skewed) : amp=158547±504, mu=92.802±0.014, sigma=0.930±0.012, skewness=-1.15±0.05
1/1 Gaussian (skewed) : amp=79559±321, mu=92.833±0.018, sigma=0.617±0.012, skewness=-0.93±0.07

Fig. 3: Power Spectral Density

path = '/content/drive/MyDrive/CityMagData/3960Hz/'
berkeley3960_day   = nuri.get_data('2016-3-14-9','2016-3-14-11',path,station=2,scalar=True,sensor='biomed',rate=3960)
berkeley3960_night = nuri.get_data('2016-3-14-2','2016-3-14-4' ,path,station=2,scalar=True,sensor='biomed',rate=3960)
brooklyn3960_day   = nuri.get_data('2018-5-7-9' ,'2018-5-7-11' ,path,station=4,scalar=True,sensor='biomed',rate=3960)
brooklyn3960_night = nuri.get_data('2018-5-7-2' ,'2018-5-7-4'  ,path,station=4,scalar=True,sensor='biomed',rate=3960)

Low and high PSDs

berkeley_day   = nuri.data_chunk(berkeley.value,[9,11])
berkeley_night = nuri.data_chunk(berkeley.value,[2,4])
brooklyn_day   = nuri.data_chunk(brooklyn.value,[9,11])
brooklyn_night = nuri.data_chunk(brooklyn.value,[2,4])

dict = {'Berkeley':{'day1Hz'      : berkeley_day,
                    'night1Hz'    : berkeley_night,
                    'day3960Hz'   : berkeley3960_day.value,
                    'night3960Hz' : berkeley3960_night.value},
        'Brooklyn':{'day1Hz'      : brooklyn_day,
                    'night1Hz'    : brooklyn_night,
                    'day3960Hz'   : brooklyn3960_day.value,
                    'night3960Hz' : brooklyn3960_night.value}
        }

import numpy
from scipy import signal
import matplotlib.pyplot as plt
fig,ax = plt.subplots(2,2,figsize=(12,6),dpi=300,sharex='row',sharey='row')
for i,name in enumerate(dict.keys()):
  for n,data in enumerate([dict[name]['day1Hz'],dict[name]['night1Hz'],dict[name]['day3960Hz'],dict[name]['night3960Hz']]):
    color = 'navy' if n in [0,2] else 'orange'
    label = 'Daytime' if n==2 else 'Nighttime' if n==3 else None
    f, Pxx_den = signal.welch(data, 1 if n<2 else 3960, nperseg=5e4 if n<2 else 1e6)
    idx = numpy.where(f>=0.1) if n in [2,3] else numpy.where(f<0.1)
    f, Pxx_den = f[idx], Pxx_den[idx]
    ax[0][i].loglog(f, Pxx_den,color=color,label=label,lw=0.5)
    ax[0][i].fill_between(f, Pxx_den,color=color,alpha=0.2)
    if n in [2,3]:
      f, Pxx_den = signal.welch(data, 3960, nperseg=1e4)
      ax[1][i].semilogy(f,Pxx_den,color=color,lw=(2 if n==2 else 1),label=label)
  # General subplot fix
  ax[0][i].axvline(0.1,color='black',lw=2,ls='dashed',alpha=0.8)
  ax[0][i].legend(loc='best',handlelength=1,handleheight=0.07)
  ax[0][i].set_xlim([1/3960,3960/2])
  ax[0][i].xaxis.tick_top()
  for n in range(1,33):
    ax[1][i].axvline(60*n,lw=0.3,ls='dashed',color='grey')
  ax[1][i].legend(loc='best',handlelength=1,handleheight=0.07)
  ax[1][i].set_xlim([0,3960/2])
  ax[1][i].set_xlabel('Frequency [Hz]')
  ax[1][i].set_yscale('log')
# General figure fix
ax[0][0].set_ylabel(r'PSD [$\mathrm{\mu T}^2/\mathrm{Hz}$]')
ax[1][0].set_ylabel(r'PSD [$\mathrm{\mu T}^2/\mathrm{Hz}$]')
plt.tight_layout()
plt.savefig('psd.png')
plt.show()

High-frequency PSD

berkeley_day   = nuri.data_chunk(berkeley.value,[5,24])
berkeley_night = nuri.data_chunk(berkeley.value,[0,5])
brooklyn_day   = nuri.data_chunk(brooklyn.value,[5,24])
brooklyn_night = nuri.data_chunk(brooklyn.value,[0,5])

dict = {'Berkeley':{'day1Hz'      : berkeley_day,
                    'night1Hz'    : berkeley_night,
                    'day3960Hz'   : berkeley3960_day.value,
                    'night3960Hz' : berkeley3960_night.value},
        'Brooklyn':{'day1Hz'      : brooklyn_day,
                    'night1Hz'    : brooklyn_night,
                    'day3960Hz'   : brooklyn3960_day.value,
                    'night3960Hz' : brooklyn3960_night.value}
        }

import numpy
from scipy import signal
import matplotlib.pyplot as plt
fig,ax = plt.subplots(2,2,figsize=(12,6),dpi=300,sharex='row',sharey='row')
for i,name in enumerate(dict.keys()):
  for n,data in enumerate([dict[name]['day3960Hz'],dict[name]['night3960Hz']]):
    color = 'navy' if n==0 else 'orange'
    label = 'Daytime' if n==0 else 'Nighttime'
    f, Pxx_den = signal.welch(data, 3960, nperseg=1e6)
    ax[0][i].loglog(f, Pxx_den,color=color,label=label,lw=0.5)
    ax[0][i].fill_between(f, Pxx_den,color=color,alpha=0.2)
    f, Pxx_den = signal.welch(data, 3960, nperseg=1e4)
    ax[1][i].semilogy(f,Pxx_den,color=color,lw=(2 if n==0 else 1),label=label)
  # General subplot fix
  ax[0][i].legend(loc='best',handlelength=1,handleheight=0.07)
  ax[0][i].set_xlim([0.1,3960/2])
  ax[0][i].set_title('Berkeley' if i==0 else 'Brooklyn',fontsize=14)
  # ax[0][i].xaxis.tick_top()
  for n in range(1,33):
    ax[1][i].axvline(60*n,lw=0.3,ls='dashed',color='grey')
  ax[1][i].legend(loc='best',handlelength=1,handleheight=0.07)
  ax[1][i].set_xlim([0,3960/2])
  ax[1][i].set_xlabel('Frequency [Hz]')
  ax[1][i].set_yscale('log')
# General figure fix
ax[0][0].set_ylabel(r'PSD [$\mathrm{\mu T}^2/\mathrm{Hz}$]')
ax[1][0].set_ylabel(r'PSD [$\mathrm{\mu T}^2/\mathrm{Hz}$]')
plt.tight_layout(h_pad=0.1)
plt.savefig('psd_high_freq.png')
plt.show()

Fig. 4: Periodic event

berkeley_day   = nuri.data_chunk(berkeley.value,[5,24])
berkeley_night = nuri.data_chunk(berkeley.value,[0,5])
brooklyn_day   = nuri.data_chunk(brooklyn.value,[5,24])
brooklyn_night = nuri.data_chunk(brooklyn.value,[0,5])

dict = {'Berkeley':{'day1Hz'      : berkeley_day,
                    'night1Hz'    : berkeley_night},
        'Brooklyn':{'day1Hz'      : brooklyn_day,
                    'night1Hz'    : brooklyn_night}
        }

def get_avg(data):
  avg_data, raw_data = [], []
  i=0
  length = 100
  while i+60*length<len(data):
    ts = data[i:i+60*length]
    b, a = signal.butter(4, [0.001/0.5], btype='low', analog=False)
    ts_filt = signal.filtfilt(b, a, ts)
    avg_data.append(ts_filt)
    raw_data.append(ts)
    i+=60*length
  avg_data = numpy.mean(numpy.array(avg_data),axis=0).flatten()
  raw_data = numpy.mean(numpy.array(raw_data),axis=0).flatten()
  return avg_data, raw_data

import numpy
from scipy import signal
import matplotlib.pyplot as plt
fig,ax = plt.subplots(2,2,figsize=(12,6),dpi=300,sharex='row')
for i,name in enumerate(dict.keys()):
  for n,data in enumerate([dict[name]['day1Hz'],dict[name]['night1Hz']]):
    color = 'navy' if n==0 else 'tomato'
    label = 'Daytime' if n==0 else 'Nighttime'
    f, Pxx_den = signal.welch(data, 1, nperseg=5e4)
    ax[0][i].loglog(f, Pxx_den,color=color,label=label,lw=0.5)
    ax[0][i].fill_between(f, Pxx_den,color=color,alpha=0.2)
  if i==0:
    y1, y2 = 2e1, 1e2
  if i==1:
    y1, y2 = 1e1, 2
  ax[0][i].annotate('',ha='center',va='bottom',xytext=(8.4e-4,y2),xy=(8.4e-4,y1),arrowprops={'facecolor':'black'})
  # General subplot fix
  ax[0][i].legend(loc='best',handlelength=1,handleheight=0.07)
  ax[0][i].set_xlim([1/3960.,0.1])
  ax[0][i].set_ylim([1e-3,1e2])
  ax[0][i].set_xlabel('Frequency [Hz]')
  ax[0][i].set_title('Berkeley' if i==0 else 'Brooklyn',fontsize=14)
  avg,raw = get_avg(dict[name]['day1Hz'] if i==0 else dict[name]['night1Hz'])
  ax[1][i].plot(numpy.arange(len(raw))/60,raw,color='grey',lw=0.2)
  ax[1][i].plot(numpy.arange(len(raw))/60,avg,color='navy' if i==0 else 'tomato',lw=2)
  ax[1][i].set_ylim([49.50,49.60] if i==0 else [92.45,92.65])
  ax[1][i].set_xlabel('Time [min]')
# General figure fix
ax[0][0].set_ylabel(r'PSD [$\mathrm{\mu T}^2/\mathrm{Hz}$]')
ax[1][0].set_ylabel(r'Magnetic Field [$\mathrm{\mu T}$]')
plt.tight_layout(h_pad=0.1)
plt.savefig('psd_low_freq.png')
plt.show()

Fig. 5: Wavelet Transform

import nuri
path = '/content/drive/MyDrive/CityMagData/1Hz/'
berkeley = nuri.get_data('2016-3-14','2016-3-15',path,station=2)
brooklyn = nuri.get_data('2018-5-7' ,'2018-5-8' ,path,station=4)
path = '/content/drive/MyDrive/CityMagData/3960Hz/'
berkeley3960 = nuri.get_data('2016-3-14-9','2016-3-14-10',path,station=2,sensor='biomed',rate=3960)
brooklyn3960 = nuri.get_data('2018-5-7-9' ,'2018-5-7-10' ,path,station=4,sensor='biomed',rate=3960)

import matplotlib.pyplot as plt
from matplotlib.ticker import LogLocator
fig = plt.figure(figsize=(12,6),dpi=300)
ax1 = fig.add_axes([0.06,0.47,0.44,0.36])
ax2 = fig.add_axes([0.51,0.47,0.44,0.36])
ax3 = fig.add_axes([0.06,0.09,0.44,0.36],sharex=ax1,sharey=ax1)
ax4 = fig.add_axes([0.51,0.09,0.44,0.36],sharex=ax2,sharey=ax2)
cax = fig.add_axes([0.06,0.85,0.89,0.03])
ax = [[ax1,ax2],[ax3,ax4]]
plt.subplots_adjust(hspace=0.01,wspace=0.01)
for i,dataset in enumerate([[berkeley,berkeley3960],[brooklyn,brooklyn3960]]):
  for n,data in enumerate(dataset):
    swap   = False   if n==0 else True
    length = 3600*24 if n==0 else 3960*60*5
    rate   = 1       if n==0 else 3960
    omega0 = 6
    dj     = 0.07    if n==0 else 0.1
    xlim   = [0,24]  if n==0 else [0,5]
    fmin   = 1/3960  if n==0 else 0.1
    vmin   = 1e-7
    vmax   = 1e4
    scales,freqs,spec = nuri.transform_data(data,length,sample_rate=rate,dj=dj,omega0=omega0)
    img = nuri.wavelet_axes(ax[i][n],scales,spec,freqs,xlim,fmin=fmin,vmin=vmin,vmax=vmax,cmap='Spectral_r',swap=swap)
  ax[i][0].set_ylabel('Frequency [Hz]',fontsize=15,labelpad=35)
  plt.setp(ax[0][i].get_xticklabels(), visible=False)
  if i==0:
    ax[1][i].set_xlabel('Monday Time [hour]',fontsize=15)
  else:
    ax[1][i].set_xlabel('Monday 9:00 to 9:05 AM [minute]',fontsize=15)
cb = plt.colorbar(img,cax=cax,orientation="horizontal")
cb.set_ticks(LogLocator())
cb.mappable.set_clim(vmin,vmax)
cb.set_label('Power $|\mathrm{W}_v|^2$ $[\mu T^2/\mathrm{Hz}]$',fontsize=15,labelpad=10)
cax.xaxis.set_label_position('top')
cax.xaxis.tick_top()
plt.savefig('wavelet.pdf')
plt.show()

Fig. 6: Brooklyn Signatures

input_data = [[100,580,nuri.get_vmr_data('/content/drive/MyDrive/CityMagData/field/VMR014_elevators.bin')],
              [250,730,nuri.get_vmr_data('/content/drive/MyDrive/CityMagData/field/VMR021_sameplatform.bin')],
              [400,880,nuri.get_vmr_data('/content/drive/MyDrive/CityMagData/field/VMR014_brooklyn_bridge.bin')],
              [100,580,nuri.get_vmr_data('/content/drive/MyDrive/CityMagData/field/VMR021_street.bin')],
              [250,730,nuri.get_vmr_data('/content/drive/MyDrive/CityMagData/field/VMR014_manhattan_bridge.bin')]]

Data selected from 2018-06-11 15:16:40.139638 to 2018-06-11 15:27:52.610085
Data selected from 2018-06-12 16:16:46.228931 to 2018-06-12 16:30:14.532916
Data selected from 2018-06-12 21:48:51.763096 to 2018-06-12 22:05:18.780926
Data selected from 2018-06-12 18:23:42.460797 to 2018-06-12 18:43:15.332219
Data selected from 2018-06-11 18:30:59.261255 to 2018-06-11 18:45:38.933017

titles = ['Elevator','Subway platform','Brooklyn bridge','Street','Manhattan bridge']

import numpy
import matplotlib.pyplot as plt
from scipy import signal
plt.style.use('seaborn')
fig,ax = plt.subplots(3,5,figsize=(12,6.5),dpi=300)
for i,(imin,imax,data) in enumerate(input_data):
  (imin,imax) = imin*200,imax*200
  data_y = data[imin:imax]
  y = data_y.value
  ax[0][i].plot(numpy.arange(imax-imin)/(200*60),y,color='black',lw=0.3)
  ax[0][i].fill_between(numpy.arange(imax-imin)/(200*60),y,color='navy',alpha=0.2)
  ax[0][i].set_xlim(0,8)
  ax[0][i].set_ylim(min(y),max(y))
  if i==0:
    ax[0][i].set_ylabel(r'Magnetic Field [$\mathrm{\mu T}$]')
  ax[0][i].set_xlabel('Time [min]')
  ax[0][i].set_title('(%i) %s' % (i+1,titles[i]))
  # Plot wavelet transform on second row
  scales,freqs,spec = nuri.transform_data(data_y,200*60*8,sample_rate=200,dj=0.1,omega0=6)
  img = nuri.wavelet_axes(ax[1][i],scales,spec,freqs,[0,8],fmin=0.005,vmin=1e-6,vmax=1e5,cmap='Spectral_r',hidey=True if i>0 else False)
  if i==0:
    ax[1][i].set_ylabel('Frequency [Hz]',labelpad=30)
  ax[1][i].set_xlabel('Time [min]')
  # Plot PSD on third row
  f, Pxx_den = signal.welch(data.value[imin:imax], 200, nperseg=500)
  ax[2][i].semilogy(f,Pxx_den,color='black',lw=0.5)
  ax[2][i].fill_between(f,Pxx_den,color='gold',alpha=0.2)
  ax[2][i].set_xlim(0,100)
  ax[2][i].set_ylim(1e-8,100)
  ax[2][i].set_xlabel('Frequency [Hz]')
  if i==0:
    ax[2][i].set_ylabel(r'PSD [$\mathrm{\mu T}^2/\mathrm{Hz}$]')
  else:
    plt.setp(ax[2][i].get_yticklabels(), visible=False)
plt.tight_layout(w_pad=0.2, h_pad=0.1)
plt.savefig('samples.pdf')
plt.show()

Fig. 7: Cross-correlation

path     = '/content/drive/MyDrive/CityMagData/'
street   = [nuri.get_vmr_data(path+'field/VMR021_street.bin',offset=-4),nuri.get_vmr_data(path+'field/VMR014_07m_street.bin',offset=-4)]
street2  = [nuri.get_vmr_data(path+'field/VMR021_street.bin',offset=-4),nuri.get_data('2018-6-12-14','2018-6-12-15',path+'3960Hz',station=4,rate=200,scalar=True,sensor='biomed')]
subway   = [nuri.get_vmr_data(path+'field/VMR014_sameplatform.bin',offset=-4),nuri.get_vmr_data(path+'field/VMR021_sameplatform.bin',offset=-4)]
building = [nuri.get_data('2017-10-9-6','2017-10-9-16',path+'1Hz',station=3,scalar=True),nuri.get_data('2018-5-14-6','2018-5-14-16',path+'1Hz',station=4,scalar=True)]

Data selected from 2018-06-12 14:23:42.460797 to 2018-06-12 14:43:15.332219
Data selected from 2018-06-12 14:32:34.238076 to 2018-06-12 14:37:27.443045
Data selected from 2018-06-12 14:23:42.460797 to 2018-06-12 14:43:15.332219
Data selected from 2018-06-12 12:19:22.026389 to 2018-06-12 12:31:27.127515
Data selected from 2018-06-12 12:16:46.228931 to 2018-06-12 12:30:14.532916

def normalize(data):
  data -= numpy.median(data)
  data /= max(abs(data))
  return data

data_list = [[200,street,0,[0.1],'low',4,2,'min',False],
             [200,street2,0,[10],'low',10,8,'min',False],
             [200,subway,0,[10],'low',10.8,2,'min',False],
             [1,building,0,[0.001],'low',60*10,60*9,'hour',True]]
titles = ['Same sidewalk','Different heights','Opposite subway platforms','Across street']

import numpy, math, datetime
from scipy import signal
import matplotlib.pyplot as plt
colors=['navy','salmon']
plt.style.use('seaborn')
fig,ax = plt.subplots(3,4,figsize=(12,6),dpi=300,sharey='row')
for i,(fs,(data1,data2),tmin,cutoff,btype,tmax,max_corr,scale,rescale) in enumerate(data_list):
  time_scale = 1 if scale=='min' else 60
  tmin = tmin if rescale else max(data1.t0.value+tmin,data2.t0.value+tmin)
  filt_data = []
  for n,data in enumerate([data1,data2]):
    imin = tmin if rescale else math.ceil((tmin-data.t0.value)/data.dt.value)
    imax = int(imin+fs*60*tmax)
    data = data.value.flatten()[imin:imax]
    x = numpy.linspace(0,round(len(data)/(60*fs*time_scale)),len(data))
    ax[0][i].plot(x,data,lw=1,color=colors[n],alpha=0.9)
    ax[0][i].set_xlim(x[0],x[-1])
    ax[0][i].set_xlabel('Time [%s]'%scale)
    ax[0][i].set_title('(%i) %s' % (i+1,titles[i]))
    # ax[0][i].set_ylim(-1,1)
    b, a = signal.butter(4, list(numpy.array(cutoff)/(0.5*fs)), btype=btype, analog=False)
    ts_filt = nuri.normalize(signal.filtfilt(b, a, data))
    ax[1][i].plot(x,ts_filt,lw=0.1,color=colors[n],alpha=0.9)
    ax[1][i].fill_between(x,ts_filt,lw=0,color=colors[n],alpha=0.7)
    ax[1][i].set_ylim(-1,1)
    ax[1][i].set_xlim(x[0],x[-1])
    ax[1][i].set_xlabel('Time [%s]'%scale)
    filt_data.append(ts_filt)
  min_len = min([len(filt) for filt in filt_data])
  lags, xcorr, _, _ = ax[2][i].xcorr(filt_data[0][:min_len],filt_data[1][:min_len],maxlags=fs*60*max_corr)
  lags = lags/(60*fs*time_scale)
  ax[2][i].clear()
  ax[2][i].fill_between(lags,xcorr,color='black')
  ax[2][i].set_ylim(-1,1)
  ax[2][i].set_xlim(lags[0],lags[-1])
  ax[2][i].set_xlabel('Lags [%s]'%scale)
ax[0][0].set_ylabel(r'Magnetic Field [$\mathrm{\mu T}$]')
ax[1][0].set_ylabel('Median-Normalized\nFiltered Data')
ax[2][0].set_ylabel('Cross Correlation')
plt.tight_layout(w_pad=0.2, h_pad=0.1)
plt.savefig('xcorr.pdf')
plt.show()