Step 1. Generate Vs30 GeoTIFF file

Download the latest Vs30 data (2023 version) from this link.

Extract files and start QGIS.

Open Project from the menu and select vs30.qgs file

There are a number of layers you can choose from.  As per Kevin Foster's paper, Geology Vs30 and Terrain Vs30 are combined into Combined Vs30. 

Let's export this to a GeoTiff. Right-click on the layer name, and Export > Save As 

Keep Raw data as Output mode, choose GeoTIFF format, and enter a file name. (Click ... button and specify the directory too). CRS is important.  For longitude and latitude, EOSGL4326 - WGS84. Other files can be left unchanged. Click OK to proceed.

Now you have obtained combined_vs30.tif 

Step 2. Extract Vs30 for given locations

First, you need a CSV file containing list of stations in "station name, longitude, latitude" format (sample_stations.vs30).

Station,lon,lat
ADCS,171.747604,-43.902401
AKSS,172.963501,-43.810902
AMBC,172.730896,-43.154701
APPS,171.567703,-42.948898
ARKS,174.944107,-41.2421
ASHS,172.595901,-43.274399
BFZ,176.246201,-40.6796

...

Use this Python code (extract_vs30.py) :  Place the GeoTIFF file and station csv file in the same directory as this code and run.

import rasterio
import numpy as np
from pathlib import Path

cwd = Path.cwd()
fname_dataset = cwd / "combined_vs30.tif"


fname_ll = "sample_stations.csv"
fname_ll = cwd / Path(fname_ll)
fname_vs30=fname_ll.with_suffix('.vs30')

ll = np.genfromtxt(fname_ll, delimiter=',', dtype=None)

n_sta = len(ll)

dataset = rasterio.open(fname_dataset)
band1 = dataset.read(1)

lat_sta = []
lon_sta = []
code_sta = []
vs30_sta = []

for ind, sta in enumerate(ll[1:]): #skip the header line
    lon = float(sta[1])
    lat = float(sta[2])
    lon_sta.append(lon)
    lat_sta.append(lat)
    code_sta.append(sta[0])
    py, px = dataset.index(lon, lat)
    vs30_sta.append(band1[py, px])

# output vs30 file
vs30file = open(fname_vs30, 'w')
for code, vs30 in zip(code_sta, vs30_sta):
    vs30file.writelines(code.decode() + ',' + str(vs30) + '\n')
vs30file.close()

It will generate sample_stations.vs30 that looks like

ADCS,302.3162
AKSS,624.8928
AMBC,224.83485
APPS,304.1254
ARKS,353.8225
ASHS,247.1665
BFZ,713.0071
BMTS,363.64764
BOWS,316.6727
BSWZ,735.06323
BWHS,205.3757
BWRS,779.0794
CACS,386.50256

...

WNAS,229.21785
WNHS,605.6527
WNKS,351.3209
WPWS,334.91064
WRCS,324.52124
WSTS,-32767.0
WTMC,328.34332
WVAS,474.6339
WVFS,324.52124
WVZ,486.30835