<jupyter_start><jupyter_text>Practical 1: Processing ground-based gravity dataIn this practical, we'll process a set of synthetic gravity surveys near the city of Cape Town, South Africa. The surveys were conducted with a LaCoste & Romberg Model G gravimeter and tied to a single base station (station 0). The synthetic data were generated based on public domain data for Southern Africa distributed by the [NOAA National Centers for Environmental Information](https://www.ngdc.noaa.gov/). Goals* Process the survey data from raw readings into absolute gravity* Merge them into a single database Data For this practical, you'll need:* A zip archive with the survey data: [`cape-town-surveys.zip`](https://github.com/leouieda/gravity-processing/raw/main/data/cape-town-surveys.zip)* The gravimeter calibration table to convert readings to mGal: [`gravimeter-scaling.csv`](https://raw.githubusercontent.com/leouieda/gravity-processing/main/data/gravimeter-scaling.csv)Download both files, unzip the archive, and place everything into the **same folder as this notebook**. Your should folder structure should look like this:```gravity-practicals/├── cape-town-surveys/│   ├── cape-town-gravity-day-1.csv│   ├── cape-town-gravity-day-2.csv│   ├── cape-town-gravity-day-3.csv│   ├── cape-town-gravity-day-4.csv│   └── cape-town-gravity-day-5.csv├── gravimeter-scaling.csv└── practical1.ipynb``` Survey dataThe survey files are in Comma Separated Values (CSV) format, which is a table of values with rows representing individual survey points and columns representing the position, time, and meter reading. For example, these are the first few lines of one of these files:```station_id,longitude,latitude,easting,northing,elevation,time_minutes,reading0,18.34444,-34.12971,255105.43,6220276.33,32.2,0.0,2555.081,18.37418,-34.19583,258037.64,6213013.01,18.4,8.0,2565.122,18.40388,-34.23972,260899.25,6208214.71,25.0,14.0,2569.493,18.41112,-34.16444,261353.99,6216582.03,228.7,23.0,2514.59```* Station ID is a unique identifier (across all surveys) for each survey point. * Latitude and longitude are geographic coordinates in degress referenced to the WGS84 ellipsoid.* Easting and northing are UTM coordinates in meters.* Elevation is the height of gravimeter above the WGS84 ellipsoid in meters.* Time are the minutes since the first reading at the base station.* Reading is the gravimeter reading (unscaled) at that stattion.The surveys were conducted in a loop, beginning and ending at the base station. Gravimeter calibration tableThe gravimeter calibration table is particular to the gravimeter used for this survey. The table was provided by the manufacturer. This is what the first few lines of the table look like:```counter_reading,value_mgal,interval_factor0000,0000.00,1.006360100,0100.64,1.006210200,0201.26,1.006090300,0301.87,1.005970400,0402.46,1.005880500,0503.05,1.005790600,0603.63,1.005700700,0704.20,1.00563```Convert the reading to mGal using the table values and the following equation:```reading_mgal = (reading - counter_reading) * interval_factor + value_mgal``` Base stationThe base station for all surveys was the same. Here are it's properties:* Station ID: 0* Absolute gravity: 979656.12 mGal  Survey locationThe surveys were caried out near Cape Town, South Africa. Most points fall around False Bay and the Cape Peninsula, stretching a bit northwards. See the Wikipedia page on the [Geology of Cape Town](https://en.wikipedia.org/wiki/Geology_of_Cape_Town) for details on the geologic context.Here is a cool thing you can do in notebooks: embed an interactive Google Map (or any HTML page, actually)! To do this, get the HTML code to embed a map by clicking on the "share" button on Google Maps. Copy the `src` part of the code and use `IPython.display.IFrame` to insert it into the notebook.<jupyter_code>from IPython.display import IFrame
IFrame(src="https://www.google.com/maps/d/embed?mid=1FunoJz1_shZ3wxfNFpg5FH33qkOg3tie", width=800, height=600)<jupyter_output><empty_output><jupyter_text>Import the required librariesTo deal with this type of tabular data, all we need are numpy and matplotlib.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Reading CSV data with numpyTo import the CSV data into a numpy array, we'll use the `numpy.loadtxt` function. By default, it expects spaces between values instead of commas, so we have to use the `delimiter` argument to specify that it needs to use commas. We also need to skip the first line of the file that includes the header (names of each column) using the `skiprows` argument.<jupyter_code>data = np.loadtxt("cape-town-surveys/cape-town-gravity-day-1.csv", delimiter=",", skiprows=1)
data<jupyter_output><empty_output><jupyter_text>---- Your turnCreate variables called `station_id`, `easting`, `northing`, `time`, and `reading` with the respective data values (columns in the `data` array).Tip: Use array slicing (like `variable[5:10, 7:8]`) to separate each column of the array.<jupyter_code>station_id = data[:, 0]
easting = data[:, 3]
northing = data[:, 4]
time = data[:, 6]
reading = data[:, 7]<jupyter_output><empty_output><jupyter_text>This would also work (`.T` transposes an array):<jupyter_code>station_id, easting, northing, time, reading = data[:, [0, 3, 4, 6, 7]].T<jupyter_output><empty_output><jupyter_text>---- Plot the survey locations on a mapNow that we have the easting and northing coordinates of each point, we can plot them on a map to see what kind of layout this survey had.<jupyter_code>plt.figure(figsize=(14, 10))
plt.plot(easting, northing, ".k")
for i in range(station_id.size):
    plt.text(easting[i], northing[i], station_id[i])
plt.xlabel("UTM Easting [m]")
plt.ylabel("UTM Northing [m]")
plt.title("Survey 1: station locations and IDs")<jupyter_output><empty_output><jupyter_text>Drift correctionFirst thing we need to do is correct for the drift. Let's plot the first and last readings to see how much drift there was.<jupyter_code>plt.figure(figsize=(10, 7))
plt.plot(time[0], reading[0], ".k")
plt.plot(time[-1], reading[-1], ".k")
plt.title("First and last reading at the base station")
plt.xlabel("time [minutes]")
plt.ylabel("gravimeter reading")<jupyter_output><empty_output><jupyter_text>Assuming that the drift is linear, we can calculate the drift factor $\alpha = \dfrac{\Delta r}{\Delta t}$<jupyter_code>drift_factor = (reading[-1] - reading[0]) / (time[-1] - time[0])
print(drift_factor)<jupyter_output>-0.029999999999999465<jupyter_text>---- Your turnApply the drift correction to the readings ($r_{corrected} = r - \alpha\ t$) and store the results in a variable called `reading_nodrift`. Make a plot of the raw readings and the corrected readings as a function of time (reading in y-axis, time in x-axis).<jupyter_code>reading_nodrift = reading - drift_factor * time
plt.figure(figsize=(10, 7))
plt.plot(time, reading, ".k", label="raw reading")
plt.plot(time, reading_nodrift, ".r", label="drift corrected")
plt.legend()
plt.xlabel("time (minutes)")
plt.ylabel("gravimeter reading")<jupyter_output><empty_output><jupyter_text>---- Convert readings to mGalThe next step is to convert the drif-corrected readings to mGal using the conversion table for this gravimeter. First, we have to load the conversion table using numpy.<jupyter_code>conversion_table = np.loadtxt("gravimeter-scaling.csv", delimiter=",", skiprows=1)
conversion_table
# Separate the columns of the table into variables
counter_reading = conversion_table[:, 0]
value_mgal = conversion_table[:, 1]
interval_factor = conversion_table[:, 2]<jupyter_output><empty_output><jupyter_text>To illustrate how this is done, let's convert the first reading to mGal. We'll use the values in the table and the formula given above for the conversion.The trick is figuring out the position (row) in the table that has the values we want for a particular reading. The reading has to be within the counter interval for that row. For example, a reading of 2710 would use the values in row 27 of the table. The position/row can be calculated using the integer division (ignoring the remainder and decimal places) of the reading by 100. With this, we can get the exact conversion factor and value in mGal from the table.<jupyter_code>reading0 = reading_nodrift[0]
print(reading0)
position = int(reading0 // 100)
print(position)
print(counter_reading[position], value_mgal[position], interval_factor[position])
reading0_mgal = (reading0 - counter_reading[position]) * interval_factor[position] + value_mgal[position]
print(reading0_mgal)<jupyter_output>2569.1108223999995<jupyter_text>Now that we know how to convert one reading, we can repeat this for all readings using a `for` loop. We can use the `range` function to generate a list of values from 0 to a given number. These will be the position of the readings in our array.<jupyter_code># Start with an empty array and fill it one element at a time
reading_mgal = np.empty(reading_nodrift.size)
for i in range(reading_nodrift.size):
    position = int(reading_nodrift[i] // 100)
    reading_mgal[i] = (reading_nodrift[i] - counter_reading[position]) * interval_factor[position] + value_mgal[position]
print(reading_mgal)<jupyter_output>[2569.1108224 2579.4451008 2584.0191248 2529.1006784 2575.7255648
 2575.464192  2555.8712848 2553.2072928 2553.3178736 2546.9343456
 2549.0152752 2546.3311776 2544.4915152 2553.5289824 2516.8865264
 2552.8051808 2562.7071888 2567.7738    2487.04556   2569.1108224]<jupyter_text>Calculate absolute gravityNow that we have our relative readings in mGal, we can tie them to the base station and calculate absolute gravity:$$ g_{station} = (r_{station} - r_{base}) + g_{base} $$The first thing we need to do is calculate the difference between the readings and the base station (station 0).<jupyter_code>reading_difference = reading_mgal - reading_mgal[0]
print(reading_difference)
gravity_base = 979656.12
gravity = reading_difference + gravity_base
print(gravity)<jupyter_output>[979656.12      979666.4542784 979671.0283024 979616.109856
 979662.7347424 979662.4733696 979642.8804624 979640.2164704
 979640.3270512 979633.9435232 979636.0244528 979633.3403552
 979631.5006928 979640.53816   979603.895704  979639.8143584
 979649.7163664 979654.7829776 979574.0547376 979656.12     ]<jupyter_text>Now we can plot the absolute gravity values on a map using the `plt.scatter` function.<jupyter_code>plt.figure(figsize=(7, 8))
plt.scatter(easting, northing, c=gravity, s=60)
plt.colorbar(label="mGal")
plt.xlabel("UTM Easting [m]")
plt.ylabel("UTM Northing [m]")
plt.title("Absolute gravity")<jupyter_output><empty_output><jupyter_text>---- Your turnRepeat the calculations above for every survey in the zip archive. Use a `for` loop to process all surveys instead of copy/pasting code. Combine the arrays for coordinates, elevation, and absolute gravity from each survey into arrays with all data from all surveys (use `numpy.concatenate` to join multiple arrays). Make a map of the absolute gravity values of the entire dataset.Tip: Use the `glob` module to get a list of all file names in a folder.<jupyter_code>import glob
# The * is a special character that means "any character in any amount"
survey_files = glob.glob("cape-town-surveys/*.csv")
print(survey_files)<jupyter_output>['cape-town-surveys/cape-town-gravity-day-5.csv', 'cape-town-surveys/cape-town-gravity-day-1.csv', 'cape-town-surveys/cape-town-gravity-day-2.csv', 'cape-town-surveys/cape-town-gravity-day-3.csv', 'cape-town-surveys/cape-town-gravity-day-4.csv']<jupyter_text>Start with empty lists for the coordinates, readings, etc. We'll store the arrays for each survey in these lists. We can then join them together using `numpy.concatenate`.<jupyter_code>latitudes = []
longitudes = []
eastings = []
northings = []
elevations = []
gravities = []
for survey_file in survey_files:
    print(f"Processing file '{survey_file}'...", end=" ")
    
    # Load the data
    data = np.loadtxt(survey_file, delimiter=",", skiprows=1)
        
    # Get the time and reading to convert them
    time = data[:, 6]
    reading = data[:, 7]

    # Drift correction
    drift_factor = (reading[-1] - reading[0]) / (time[-1] - time[0])
    reading_nodrift = reading - drift_factor * time

    # Convert to mGal
    reading_mgal = np.empty(reading_nodrift.size)
    for i in range(reading_nodrift.size):
        position = int(reading_nodrift[i] // 100)
        reading_mgal[i] = (reading_nodrift[i] - counter_reading[position]) * interval_factor[position] + value_mgal[position]
        
    # Calculate absolute gravity
    reading_difference = reading_mgal - reading_mgal[0]
    gravity = reading_difference + gravity_base

    # Append the arrays for this survey to the lists
    longitudes.append(data[:, 1])
    latitudes.append(data[:, 2])
    eastings.append(data[:, 3])
    northings.append(data[:, 4])
    elevations.append(data[:, 5])
    gravities.append(gravity)
    
    print("Done")
    
longitude = np.concatenate(longitudes)
latitude = np.concatenate(latitudes)
easting = np.concatenate(eastings)
northing = np.concatenate(northings)
elevation = np.concatenate(elevations)
gravity = np.concatenate(gravities)
plt.figure(figsize=(7, 8))
plt.scatter(easting, northing, c=gravity, s=60)
plt.colorbar(label="mGal")
plt.xlabel("UTM Easting [m]")
plt.ylabel("UTM Northing [m]")
plt.title("Absolute gravity from all surveys")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10483271375464684, 'fraction_numerical': 0.06951672862453531, 'mean_word_length': 5.380929791271347, 'pattern_counts': {'":': 0, '<': 60, '<?xml version=': 0, '>': 60, 'https://': 5, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '218649', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 4940, 'n_tokens_neox': 4229, 'n_words': 1404, 'path': 'practicals/practical1-solution.ipynb', 'repo_name': 'leouieda/gravity-processing', 'size': '13450'}