Add Git LFS tracking for large files

Dockerfile

@@ -0,0 +1,22 @@
+FROM python:3.11
+
+# Set up a new user named "user" with user ID 1000 for permission
+RUN useradd -m -u 1000 user
+# Switch to the "user" user
+USER user
+# Set home to the user's home directory
+ENV HOME=/home/user \
+	PATH=/home/user/.local/bin:$PATH
+
+# Upgreade pip
+RUN pip install --no-cache-dir --upgrade pip
+
+COPY --chown=user requirements.txt .
+
+# Install requirements
+RUN pip install --no-cache-dir --upgrade -r requirements.txt
+
+COPY --chown=user *.py *.css  /
+COPY --chown=user helpers/* /helpers/
+
+ENTRYPOINT ["solara", "run", "app.py", "--host=0.0.0.0", "--port", "7860", "--production"]

README.md

@@ -1,10 +1,11 @@
 ---
-title: Example1
-emoji: ⚡
-colorFrom: pink
-colorTo: indigo
+title: Major TOM Viewer
+emoji: 🌍
+colorFrom: gray
+colorTo: gray
 sdk: docker
 pinned: false
+short_description: Quick View of Samples in the MajorTOM-Core Dataset
 ---
 
-Check out the configuration reference at https://huggingface.co/docs/hub/spaces-config-reference
+Check out the configuration reference at https://huggingface.co/docs/hub/spaces-config-reference

app.py

@@ -0,0 +1,92 @@
+import sys
+import os
+import leafmap
+from helpers.helpers_grid import *
+from helpers.helpers_functional import *
+from leafmap.toolbar import change_basemap
+from IPython.display import display
+import ipywidgets
+import solara
+
+instructions_top = '''
+![image/png](https://cdn-uploads.huggingface.co/production/uploads/6304c06eeb6d777a838eab63/BJKsLwX0GG4W3-gdf40TJ.png)
+# Dataset Viewer
+This app provides a way of exploring samples present in the MajorTOM-Core dataset.
+It contains nearly every piece of Earth captured by ESA Sentinel-2 satellite as well as a large fraction of paired Sentinel-1 data from a similar time period.
+### Instructions
+To find a sample, navigate on the map to a place of interest. Click `Find Sample` to find a dataset sample that contains the central pixel of your current view.
+'''
+
+instructions_bottom = '''
+<details><summary><strong>🏝 Couldn't find a sample? See this figure of global coverage:</strong></summary>
+    <img src='https://cdn-uploads.huggingface.co/production/uploads/6304c06eeb6d777a838eab63/2KTarfsM0a1dNYEbXriUH.png' />
+</details>
+'''
+
+image_data = solara.reactive(Image.new('RGB',(1068,1068)))
+center = solara.reactive((20, 0))
+gridcell = solara.reactive('')
+timestamp = solara.reactive('')
+zoom= solara.reactive(4)
+source = solara.reactive('S2-L2A')
+
+@solara.component
+def Page():
+    with solara.Column():
+
+        with solara.Card(margin=10):
+            solara.Markdown(instructions_top)
+            solara.Button(label="Website",
+                  icon_name="mdi-map-legend",
+                  attributes={"href": 'https://www.huggingface.co/Major-TOM', "target": "_blank"},
+                  text=True,
+                  outlined=True)
+    
+            solara.Button(label="arXiv Paper",
+                          icon_name="mdi-script-text",
+                          attributes={"href": 'https://www.arxiv.org/abs/2402.12095', "target": "_blank"},
+                          text=True,
+                          outlined=True)
+            solara.Markdown(instructions_bottom)
+        
+        def update_image():
+            ret = map_to_image(m, return_centre=True, return_gridcell=True, return_timestamp=True, source=source.value)
+            if ret is not None:
+                image_data.value, center.value, gridcell.value, timestamp.value = ret
+                zoom.value=12
+            else:
+                image_data.value = Image.new('RGB',(1068,1068))
+                center.value = (20,0)
+                zoom.value = 4
+
+        def update_source(val):
+            source.value = val
+            update_image()
+        
+        with solara.Card(margin=10):
+            with solara.ColumnsResponsive(default=12, small=12, medium=6, large=6, xlarge=6):
+                with solara.Column(align='center'):
+                    m = leafmap.Map(
+                        height=560,
+                        width=560,
+                        layout=ipywidgets.Layout(max_width='60vw', max_height='80vh'),
+                        zoom=zoom.value,
+                        center=center.value,
+                        draw_control=False,
+                        measure_control=False,
+                        fullscreen_control=False,
+                        toolbar_control=False,
+                        attribution_control=True,
+                    )
+                    display(m)
+                    button = solara.Button("Find Sample", on_click=update_image, height='80px')                    
+                    solara.Select(label="Processing Level", value=source, values=['S2-L2A','S2-L1C', 'S1-RTC'], on_value=update_source)
+                with solara.Column(align='center'):
+                    output = solara.Image(image_data.value)
+                    solara.Markdown('''        
+                    |                     |      |
+                    |--------------------:|:----:|
+                    | Latitude, Longitude | {:.3f}, {:.3f} |
+                    |       MajorTOM Grid | {} |
+                    | Timestamp| {} |
+                    '''.format(center.value[0], center.value[1], gridcell.value, timestamp.value))

helpers/helpers_functional.py

@@ -0,0 +1,120 @@
+from fsspec.parquet import open_parquet_file 
+import fsspec
+import pyarrow.parquet as pq
+from .helpers_grid import *
+import pandas as pd
+from io import BytesIO
+import os
+from PIL import Image
+import datetime
+
+# GLOBAL VARIABLES
+if os.path.isfile('helpers/s2l2a_metadata.parquet'):
+    l2a_meta_path = 'helpers/s2l2a_metadata.parquet'
+else:
+    DATASET_NAME = 'Major-TOM/Core-S2L2A'
+    l2a_meta_path = 'https://huggingface.co/datasets/{}/resolve/main/metadata.parquet'.format(DATASET_NAME)
+
+if os.path.isfile('helpers/s2l1c_metadata.parquet'):
+    l1c_meta_path = 'helpers/s2l1c_metadata.parquet'
+else:
+    DATASET_NAME = 'Major-TOM/Core-S2L1C'
+    l1c_meta_path = 'https://huggingface.co/datasets/{}/resolve/main/metadata.parquet'.format(DATASET_NAME)
+
+if os.path.isfile('helpers/s1rtc_metadata.parquet'):
+    rtc_meta_path = 'helpers/s1rtc_metadata.parquet'
+else:
+    DATASET_NAME = 'Major-TOM/Core-S1RTC'
+    rtc_meta_path = 'https://huggingface.co/datasets/{}/resolve/main/metadata.parquet'.format(DATASET_NAME)
+
+grid = Grid(10, latitude_range=(-90,90), longitude_range=(-180,180))
+l2a_df = pd.read_parquet(l2a_meta_path)
+l1c_df = pd.read_parquet(l1c_meta_path)
+rtc_df = pd.read_parquet(rtc_meta_path)
+
+df_dict = {
+    'S2-L2A' : l2a_df,
+    'S2-L1C' : l1c_df,
+    'S1-RTC' : rtc_df
+}
+
+def pretty_date(input):
+    template = '%Y%m%dT%H%M%S' if 'T' in input else '%Y%m%d%H%M%S'
+    return datetime.datetime.strptime(input, template).strftime('%H:%M:%S - %d %b %Y')
+
+# HELPER FUNCTIONS
+def gridcell2ints(grid_string):
+    up = int(grid_string.split('_')[0][:-1]) * (2*int(grid_string.split('_')[0][-1]=='U') - 1) # +ve if up
+    right = int(grid_string.split('_')[1][:-1]) * (2*int(grid_string.split('_')[1][-1]=='R') - 1) # +ve if R
+
+    return up, right
+
+def row2image(parquet_url, parquet_row, fullrow_read=True):
+
+    if fullrow_read:
+        # option 1
+        f=fsspec.open(parquet_url)
+        temp_path = f.open()
+    else:
+        # option 2
+        temp_path = open_parquet_file(parquet_url,columns = ["thumbnail"])
+    
+    with pq.ParquetFile(temp_path) as pf:
+        first_row_group = pf.read_row_group(parquet_row, columns=['thumbnail'])
+
+    stream = BytesIO(first_row_group['thumbnail'][0].as_py())
+    return Image.open(stream)
+
+def row2s2(parquet_url, parquet_row, s2_bands = ["B04", "B03", "B02"]):
+    with open_parquet_file(parquet_url,columns = s2_bands) as f:
+        with pq.ParquetFile(f) as pf:
+            first_row_group = pf.read_row_group(parquet_row, columns=s2_bands)
+
+    return first_row_group
+
+def cell2row(grid_string, meta_df, return_row = False):
+    row_U, col_R = gridcell2ints(grid_string)
+    R = meta_df.query('grid_row_u == {} & grid_col_r == {}'.format(row_U, col_R))
+
+    if not R.empty:
+        if return_row:
+            return R.parquet_url.item(), R.parquet_row.item(), R
+        else:
+            return R.parquet_url.item(), R.parquet_row.item()
+    else:
+        return None
+
+def map_to_image(map, return_centre=False, return_gridcell=False, return_timestamp=False, source='S2-L2A'):
+
+    try:
+        # 1. get bounds
+        bbox = map.get_bbox()
+        center = [(bbox[3]+bbox[1])/2, (bbox[2]+bbox[0])/2]
+    except:
+        return None
+
+    # 2. translate coordinate to major-tom tile
+    rows, cols = grid.latlon2rowcol([center[0]], [center[1]])
+
+    # 3. translate major-tom cell to row in parquet
+    df = df_dict[source]
+    row = cell2row("{}_{}".format(rows[0],cols[0]), df, return_row = True)
+
+    if row is not None:
+        parquet_url, parquet_row, meta_row = row
+        print(meta_row)
+        img = row2image(parquet_url, parquet_row)
+        # 4. acquire image # X. update map
+        lat, lon = meta_row.centre_lat.item(), meta_row.centre_lon.item()
+
+        ret = [img]
+        if return_centre:
+            ret.append((lat,lon))
+        if return_gridcell:
+            ret.append(meta_row.grid_cell.item())
+        if return_timestamp:
+            ret.append(pretty_date(meta_row.timestamp.item()))
+            
+        return ret
+    else:
+        return None

helpers/helpers_grid.py

@@ -0,0 +1,256 @@
+import numpy as np
+import math
+import pandas as pd
+import geopandas as gpd
+from shapely.geometry import LineString, Polygon
+from tqdm import tqdm
+
+
+
+class Grid():
+
+    RADIUS_EQUATOR = 6378.137 # km
+
+    def __init__(self,dist,latitude_range=(-85,85),longitude_range=(-180,180),utm_definition='bottomleft'):
+        self.dist = dist
+        self.latitude_range = latitude_range
+        self.longitude_range = longitude_range
+        self.utm_definition = utm_definition
+        self.rows,self.lats = self.get_rows()
+        self.points, self.points_by_row = self.get_points()
+        
+    def get_rows(self):
+
+        # Define set of latitudes to use, based on the grid distance
+        arc_pole_to_pole = math.pi * self.RADIUS_EQUATOR
+        num_divisions_in_hemisphere = math.ceil(arc_pole_to_pole / self.dist)
+
+        latitudes = np.linspace(-90, 90, num_divisions_in_hemisphere+1)[:-1]
+        latitudes = np.mod(latitudes, 180) - 90
+
+        # order should be from south to north
+        latitudes = np.sort(latitudes)
+
+        zeroth_row = np.searchsorted(latitudes,0)
+
+        # From 0U-NU and 1D-ND
+        rows = [None] * len(latitudes)
+        rows[zeroth_row:] = [f'{i}U' for i in range(len(latitudes)-zeroth_row)]
+        rows[:zeroth_row] = [f'{abs(i-zeroth_row)}D' for i in range(zeroth_row)]
+
+        # bound to range
+        idxs = (latitudes>=self.latitude_range[0]) * (latitudes<=self.latitude_range[1])
+        rows,latitudes = np.array(rows), np.array(latitudes)
+        rows,latitudes = rows[idxs],latitudes[idxs]
+
+        return rows,latitudes
+
+    def get_circumference_at_latitude(self,lat):
+
+        # Circumference of the cross-section of a sphere at a given latitude
+
+        radius_at_lat = self.RADIUS_EQUATOR * math.cos(lat * math.pi / 180)
+        circumference = 2 * math.pi * radius_at_lat
+
+        return circumference
+
+    def subdivide_circumference(self,lat,return_cols=False):
+        # Provide a list of longitudes that subdivide the circumference of the earth at a given latitude
+        # into equal parts as close as possible to dist
+
+        circumference = self.get_circumference_at_latitude(lat)
+        num_divisions = math.ceil(circumference / self.dist)
+        longitudes = np.linspace(-180,180, num_divisions+1)[:-1]
+        longitudes = np.mod(longitudes, 360) - 180
+        longitudes = np.sort(longitudes)
+
+
+        if return_cols:
+            cols = [None] * len(longitudes)
+            zeroth_idx = np.where(longitudes==0)[0][0]
+            cols[zeroth_idx:] = [f'{i}R' for i in range(len(longitudes)-zeroth_idx)]
+            cols[:zeroth_idx] = [f'{abs(i-zeroth_idx)}L' for i in range(zeroth_idx)]
+            return np.array(cols),np.array(longitudes)
+
+        return np.array(longitudes)
+
+    def get_points(self):
+        
+        r_idx = 0
+        points_by_row = [None]*len(self.rows)
+        for r,lat in zip(self.rows,self.lats):
+            point_names,grid_row_names,grid_col_names,grid_row_idx,grid_col_idx,grid_lats,grid_lons,utm_zones,epsgs = [],[],[],[],[],[],[],[],[]
+            cols,lons = self.subdivide_circumference(lat,return_cols=True)
+
+            cols,lons = self.filter_longitude(cols,lons)
+            c_idx = 0
+            for c,lon in zip(cols,lons):
+                point_names.append(f'{r}_{c}')
+                grid_row_names.append(r)
+                grid_col_names.append(c)
+                grid_row_idx.append(r_idx)
+                grid_col_idx.append(c_idx)
+                grid_lats.append(lat)
+                grid_lons.append(lon)
+                if self.utm_definition == 'bottomleft':
+                    utm_zones.append(get_utm_zone_from_latlng([lat,lon]))
+                elif self.utm_definition == 'center':
+                    center_lat = lat + (1000*self.dist/2)/111_120
+                    center_lon = lon + (1000*self.dist/2)/(111_120*math.cos(center_lat*math.pi/180))
+                    utm_zones.append(get_utm_zone_from_latlng([center_lat,center_lon]))
+                else:
+                    raise ValueError(f'Invalid utm_definition {self.utm_definition}')
+                epsgs.append(f'EPSG:{utm_zones[-1]}')
+
+                c_idx += 1
+            points_by_row[r_idx] = gpd.GeoDataFrame({
+                'name':point_names,
+                'row':grid_row_names,
+                'col':grid_col_names,
+                'row_idx':grid_row_idx,
+                'col_idx':grid_col_idx,
+                'utm_zone':utm_zones,
+                'epsg':epsgs
+            },geometry=gpd.points_from_xy(grid_lons,grid_lats))
+            r_idx += 1
+        points = gpd.GeoDataFrame(pd.concat(points_by_row))
+        # points.reset_index(inplace=True,drop=True)
+        return points, points_by_row
+
+    def group_points_by_row(self):
+        # Make list of different gdfs for each row
+        points_by_row = [None]*len(self.rows)
+        for i,row in enumerate(self.rows):
+            points_by_row[i] = self.points[self.points.row==row]
+        return points_by_row
+
+    def filter_longitude(self,cols,lons):
+        idxs = (lons>=self.longitude_range[0]) * (lons<=self.longitude_range[1])
+        cols,lons = cols[idxs],lons[idxs]
+        return cols,lons
+
+    def latlon2rowcol(self,lats,lons,return_idx=False):
+        """
+        Convert latitude and longitude to row and column number from the grid
+        """
+        # Always take bottom left corner of grid cell
+        rows = np.searchsorted(self.lats,lats)-1
+
+        # Get the possible points of the grid cells at the given latitude
+        possible_points = [self.points_by_row[row] for row in rows]
+
+        # For each point, find the rightmost point that is still to the left of the given longitude
+        cols = [poss_points.iloc[np.searchsorted(poss_points.geometry.x,lon)-1].col for poss_points,lon in zip(possible_points,lons)]
+        rows = self.rows[rows]
+
+        if return_idx:
+            # Get the table index for self.points with each row,col pair in rows, cols
+            idx = [self.points[(self.points.row==row) & (self.points.col==col)].index.values[0] for row,col in zip(rows,cols)]
+            return rows,cols,idx
+        return rows,cols
+
+    def rowcol2latlon(self,rows,cols):
+        point_geoms = [self.points.loc[(self.points.row==row) & (self.points.col==col),'geometry'].values[0] for row,col in zip(rows,cols)]
+        lats = [point.y for point in point_geoms]
+        lons = [point.x for point in point_geoms]
+        return lats,lons
+
+    def get_bounded_footprint(self,point,buffer_ratio=0):
+        # Gets the polygon footprint of the grid cell for a given point, bounded by the other grid points' cells.
+        # Grid point defined as bottom-left corner of polygon. Buffer ratio is the ratio of the grid cell's width/height to buffer by.
+
+        bottom,left = point.geometry.y,point.geometry.x
+        row = point.row
+        row_idx = point.row_idx
+        col_idx = point.col_idx
+        next_row_idx = row_idx+1
+        next_col_idx = col_idx+1
+
+        if next_row_idx >= len(self.lats): # If at top row, use difference between top and second-to-top row for height
+            height = (self.lats[row_idx] - self.lats[row_idx-1])
+            top = self.lats[row_idx] + height
+        else:
+            top = self.lats[next_row_idx]
+        
+        max_col = len(self.points_by_row[row].col_idx)-1
+        if next_col_idx > max_col: # If at rightmost column, use difference between rightmost and second-to-rightmost column for width
+            width = (self.points_by_row[row].iloc[col_idx].geometry.x - self.points_by_row[row].iloc[col_idx-1].geometry.x)
+            right = self.points_by_row[row].iloc[col_idx].geometry.x + width
+        else:
+            right = self.points_by_row[row].iloc[next_col_idx].geometry.x
+
+        # Buffer the polygon by the ratio of the grid cell's width/height
+        width = right - left
+        height = top - bottom
+
+        buffer_horizontal = width * buffer_ratio
+        buffer_vertical = height * buffer_ratio
+
+        new_left = left - buffer_horizontal
+        new_right = right + buffer_horizontal
+
+        new_bottom = bottom - buffer_vertical
+        new_top = top + buffer_vertical
+
+        bbox = Polygon([(new_left,new_bottom),(new_left,new_top),(new_right,new_top),(new_right,new_bottom)])
+
+        return bbox
+    
+
+def get_utm_zone_from_latlng(latlng):
+    """
+    Get the UTM ZONE from a latlng list.
+
+    Parameters
+    ----------
+    latlng : List[Union[int, float]]
+        The latlng list to get the UTM ZONE from.
+
+    return_epsg : bool, optional
+        Whether or not to return the EPSG code instead of the WKT, by default False
+
+    Returns
+    -------
+    str
+        The WKT or EPSG code.
+    """
+    assert isinstance(latlng, (list, np.ndarray)), "latlng must be in the form of a list."
+
+    zone = math.floor(((latlng[1] + 180) / 6) + 1)
+    n_or_s = "S" if latlng[0] < 0 else "N"
+
+    false_northing = "10000000" if n_or_s == "S" else "0"
+    central_meridian = str(zone * 6 - 183)
+    epsg = f"32{'7' if n_or_s == 'S' else '6'}{str(zone)}"
+
+    return epsg
+
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+
+    dist = 100
+    grid = Grid(dist,latitude_range=(10,70),longitude_range=(-30,60))
+
+    from pprint import pprint
+
+    test_lons = np.random.uniform(-20,50,size=(1000))
+    test_lats = np.random.uniform(12,68,size=(1000))
+
+    test_rows,test_cols = grid.latlon2rowcol(test_lats,test_lons)
+    test_lats2,test_lons2 = grid.rowcol2latlon(test_rows,test_cols)
+
+    print(test_lons[:10])
+    print(test_lats[:10])
+    print(test_rows[:10])
+    print(test_cols[:10])
+
+    # Make line segments from the points to their corresponding grid points
+    lines = []
+    for i in range(len(test_lats)):
+        lines.append([(test_lons[i],test_lats[i]),(test_lons2[i],test_lats2[i])])
+
+    lines = gpd.GeoDataFrame(geometry=gpd.GeoSeries([LineString(line) for line in lines])) 
+
+    lines.to_file(f'testlines_{dist}km.geojson',driver='GeoJSON')
+    grid.points.to_file(f'testgrid_{dist}km.geojson',driver='GeoJSON')

helpers/s1rtc_metadata.parquet

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:013fe809ebf0e149d339e7e4b265e88787858092012eb8a422e6c3a1f546c117
+size 79499536

helpers/s2l1c_metadata.parquet

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e19d451d89510923c4bc3d8acf8c45985903a73c89539e431292d6226a4b5ddc
+size 171721623

helpers/s2l2a_metadata.parquet

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9a96c80bc43cb841b8400b05e80f4b477453b51a1e6833821333e3c11831e78b
+size 173048695

requirements.txt

@@ -0,0 +1,13 @@
+leafmap
+mapwidget
+solara
+geopandas
+fsspec
+pyarrow
+shapely
+pandas
+numpy
+tqdm
+requests
+aiohttp
+pydantic<2.0