Extract values for a given latitude/longitude

[quinten]

You can use some simple Python code to extract values from an ACOLITE NetCDF output.

First import the acolite code (cloned from https://github.com/acolite/acolite). In this example acolite is cloned into the git/acolite directory in the user's home folder:

Code: Select all

## add acolite to path & import
import sys, os, glob
user_home = os.path.expanduser("~")
sys.path.append('{}/git/acolite'.format(user_home))
import acolite as ac
import numpy as np

Then specify the path to the NetCDF file output by ACOLITE and the lat/lon coordinates of the point:

Code: Select all

## ACOLITE NetCDF file
ncf = '/path/to/output.nc'

## station lat/lon
st_lat = 
st_lon = 

Read the file and get the data values for the given location:

Code: Select all

## get datasets present in file
datasets = ac.shared.nc_datasets(ncf)

## read lat/lon
lon = ac.shared.nc_data(ncf, 'lon')
lat = ac.shared.nc_data(ncf, 'lat')

## find location in pixel coordinates
londiff = np.abs(lon-st_lon)    
latdiff = np.abs(lat-st_lat)
diff = np.power(np.power(londiff,2)+np.power(latdiff,2),0.5)
x,y = np.where(diff == diff.min())
x = x[0]
y = y[0]
print(st_lon, st_lat)
print(lon[x,y], lat[x,y])

## read data
extraction = {ds: ac.shared.nc_data(ncf, ds, sub=[y,x,1,1])[0][0] for ds in datasets}

## print results
for ds in extraction:
    print(ds, extraction[ds])

Note that this code does not check whether your lat/lon are actually within the image, but it should get you started with simple data extraction. You can also use the sub parameter to extract a box of values, e.g. 3x3 pixels would be sub=[y-1, x-1, 3, 3]

Quinten

[IsabelBrand]

Hi Quinten,

Great post!

Could you provide a version of this code where we could input a vector (.shp) with a few coordinates and to loop through a bunch of .nc?

Thanks in advance,


Isabel

[quinten]

Hi Isabel

If you can send me an example of the shape file I can take a look. I am not sure how these are organised.

Quinten

[quinten]

I have added a function to the ACOLITE GitHub code for extracting lat/lon points from the ACOLITE NetCDF outputs. Make sure to pull in changes from GitHub and change the paths in the example:

Code: Select all

## add acolite to path & import
import sys, os
sys.path.append('{}/git/acolite'.format(os.path.expanduser("~")))
import acolite as ac

## ACOLITE NetCDF file
ncf = '/path/to/output.nc'

## station lat/lon
st_lat = 
st_lon = 

## extract data from pixel containing station
ext1x1 = ac.shared.nc_extract_point(ncf, st_lon, st_lat, extract_datasets = None, shift_edge = False, box_size = 1)

## extract data from 5x5 box around station
ext5x5 = ac.shared.nc_extract_point(ncf, st_lon, st_lat, extract_datasets = None, shift_edge = False, box_size = 5)

The function returns a dict with "data" and also includes "mean", "median", "std", "n" when a box is extracted. "gatts" are the NetCDF global attributes.

Code: Select all

print(ext1x1.keys())

gives:

Code: Select all

dict_keys(['gatts', 'data'])

and

Code: Select all

print(ext5x5.keys()) gives

gives:

Code: Select all

dict_keys(['gatts', 'data', 'mean', 'std', 'median', 'n'])

Quinten

[IsabelBrand]

Hi Quinten,

I was trying the above instructions to extract SPM data for a couple of points, but I got the error below.
------------------------------------------------------------------------------------------------------------------------
AttributeError: module 'acolite.shared' has no attribute 'nc_extract_point'
------------------------------------------------------------------------------------------------------------------------
I only got outputs after adding the original code from https://github.com/acolite/acolite/blob ... t_point.py
after import acolite as ac
------------------------------------------------------------------------------------------------------------------------
## add acolite to path & import
import sys, os
sys.path.append('{}/git/acolite'.format(os.path.expanduser("~")))
import acolite as ac
------------------------------------------------------------------------------------------------------------------------
It worked, but I believe is not the ideal solution.

My testing scene: S2A_MSIL1C_20160501T105310_N0201_R051_T31UET_20160501T194451.SAFE
and points:
sites,lon,lat
1 4.287550775 51.89861203
2 4.267974031 51.90524905
3 4.248327738 51.91251089
4 4.223552856 51.92841046
5 4.185210008 51.94353363
6 3.981409396 52.00930478

Any suggestions on this error is much appreciated.

Cheers,

Isabel

[quinten]

Hi Isabel

Did you update the ACOLITE code base since this script was added in March or so?

Quinten

[IsabelBrand]

Hi Quinten,

Thanks! I got no errors after updating ACOLITE; however, I do have more questions on this matter.

1) I found it rather complicated to extract only ['mean', 'std', 'median', 'n'] from the dictionary. Is it nested?
2) How does it work if one wants to extract values for more than one location from multiple files (.nc)?

Best regards,

Isabel

[quinten]

Hi Isabel

The return value is a Python dictionary with the variables stored under the different statistics (e.g. mean and std), but you can reorder this to store all statistics per variable:

Code: Select all

## extract st_lon, st_lat from ncf
ext = ac.shared.nc_extract_point(ncf, st_lon, st_lat, extract_datasets = None, shift_edge = False, box_size = 5)
## reorder extracted statistics per dataset
ext_per_ds = {ds: {'mean':ext['mean'][ds], 'std':ext['std'][ds], 'median':ext['median'][ds],'n':ext['n'][ds]} for ds in ext['datasets']}

Printing a single key then gives for example, print(ext_per_ds['rhot_492']):

Code: Select all

{'mean': 0.29368597, 'std': 0.0022542998, 'median': 0.2937, 'n': 121}

If you want to do multiple stations and files you'll need to construct a nested loop (1) over files and (2) over stations. For example, set up variables:

Code: Select all

## box dimensions in pixels
box = 11

## ACOLITE file type to extract data from L1R L2R L2W ST
file_type = 'L2W'

## directory of files
idir = '/path/to/your/ACOLITE/output/directory/'

## set up list of stations
stations = [{'name': 'station1', 'st_lon':2.9074, 'st_lat':51.2555},
                {'name': 'station2', 'st_lon':2.9340, 'st_lat':51.2494,},
                {'name': 'station3', 'st_lon':2.9198, 'st_lat':51.2395}]

and extract the data:

Code: Select all

## Import ACOLITE code
import sys, os
user_home = os.path.expanduser("~")
sys.path.append(user_home+'/git/acolite')
import acolite as ac

## search files files
import glob
files = glob.glob('{}/*_{}.nc'.format(idir, file_type))
files.sort()

## run through files and stations
results = {}
for file in files:
    print('Extracting {}'.format(file))
    bn = os.path.basename(file)
    results[bn] = {}
    
    for st in stations:
        print('Extracting {}'.format(st['name']))
        ext = ac.shared.nc_extract_point(file, st['st_lon'], st['st_lat'], 
                                         extract_datasets = None, shift_edge = False, box_size = box)
        results[bn][st['name']] = {ds: {'mean':ext['mean'][ds], 
                                        'std':ext['std'][ds],
                                        'median':ext['median'][ds],
                                        'n':ext['n'][ds]} for ds in ext['mean']}

You could also replace the station dictionary with lists you index and files by a list of files:

Code: Select all

## Import ACOLITE code
import sys, os
user_home = os.path.expanduser("~")
sys.path.append(user_home+'/git/acolite')
import acolite as ac

files = ['/path/to/file1', '/path/to/file2', '/path/to/file3']
stations = ['station1', 'station2', 'station3']
lats = [51.2555, 51.2494, 51.2395]
lons = [2.9074, 2.9340, 2.9198]

results = {}
for file in files:
    print('Extracting {}'.format(file))
    bn = os.path.basename(file)
    results[bn] = {}
    
    for si, st in enumerate(stations):
        print('Extracting {}'.format(st))
        ext = ac.shared.nc_extract_point(file, lons[si], lats[si], 
                                         extract_datasets = None, shift_edge = False, box_size = box)
        results[bn][st] = {ds: {'mean':ext['mean'][ds], 
                                        'std':ext['std'][ds],
                                        'median':ext['median'][ds],
                                        'n':ext['n'][ds]} for ds in ext['mean']}

I hope this helps!

Quinten

[IsabelBrand]

Hi Quinten,

Thank you for your support. It works perfectly.

Cheers,

Isabel