MERIS TSM western US coast

[lopezama]

Hello! I would like to download the full resolution MERIS Level 2 total suspended matter (TSM) data product for the entire western US coast (shoreline to 5km offshore) for the full 2002-2012 data record. I have tried using the following data portals, but none of these actually allow me to download the specific data I am targeting. Does anyone know of a data download platform that more easily allows querying Full Resolution Level 2 MERIS TSM data by ROI for the entire data record?

Envisat MERIS https://meris-ds.eo.esa.int/oads/access/
Allows you to roughly define ROI but for some reason the interface doesn't allow the data acquisition date to be adjusted beyond the most recent week or so (ex. right now it only allows dates in Jan 2025 to be selected).

EarthData Search https://search.earthdata.nasa.gov/search
Allows you to download all Level 2 full resolution MERIS data for full data record. Cannot isolate just the TSM product, cannot define ROI then download data by ROI (must download full swath).

Ocean Color Web interface https://oceancolor.gsfc.nasa.gov/cgi/browse.pl?sen=amod
Allows you to roughly define ROI and access full data record, but does not include the full resolution MERIS TSM product.

Any help is much appreciated. Thanks!
Mandy
amanda.m.lopez@jpl.nasa.gov

[OB.DAAC - SeanBailey]

We cannot provide support for the ESA data distribution site, for that I suggest contacting the ESA services desk.

The L2 data distributed by the Ocean Color Web interface does not include a TSM product as it is not the ESA product, but rather an Ocean Biology Processing Group (OBPG) produced product that is consistent with data produced for the other sensors supported by the OBPG.

The ESA-produced data distributed by LAADS (and searchable via Earthdata Search) do not have the same level of service as other NASA-produced because the SAFE format is a zipped directory of individual netCDF files and an extraction service has not been implemented for it.

Regards,
Sean

[lopezama]

Hi Sean,
Thank you for the prompt response. Just a couple follow up questions-
1) Ocean Color Web only provides the ESA data products and the TSM product is a non-ESA product generated by the OBPG. Therefore, TSM will not be made available via Ocean Color Web. Is my understanding correct?

2) For the Level 2 Full Resolution MERIS data on Earthdata Search, is there a way to specify an ROI and only download data for that ROI? Even if there is no way to specifically download the TSM data, it would be easier to download smaller files versus the large swaths.

Thank you,
Mandy

[OB.DAAC - SeanBailey]

Mandy,

1) The Ocean Color Web only provides the OBPG-produced L2 data products for MERIS and the TSM product is an ESA product NOT generated by the OBPG. The OB.DAAC doesn't have an equivalent TSM product for any of our supported missions.

2) No, I do not believe there is any extraction capability for the ESA-produced L2 files distributed by LAADS (and available via Earthdata Search) - spatial or otherwise.

Just FYI, the OBPG-produced L2s are also available on Earthdata Search, but as per #1, don't include a TSM product.

Sorry,

Sean

[lopezama]

Thanks for your assistance. I moved forward with using the MERIS Level 2 “Full Resolution Full Swath Geophysical Product for Ocean, Land and Atmosphere” data from EarthData Search which has a total suspended matter product included. I'm struggling with georeferencing the tsm_nn.nc file using the geo_coordinates.nc file. I've been able to georeference and convert to geotiff but the final result plots in a slightly incorrect location. I have searched for resources to help with this (ex. https://help.marine.copernicus.eu/en/articles/5029956-how-to-convert-netcdf-to-geotiff), but have not found something that works well.
Any suggestions?

Mandy

[OB.DAAC - SeanBailey]

Mandy,
Did you map the data to a geotiff supported projection? The L2 netCDF files, though geolocated, are not projected.

Sean

[lopezama]

Hi Sean,
Below is my python code for an example file downloaded from EarthData Search EN1_MDSI_MER_FRS_2P_20120402T180502_20120402T180811_052785_0185_20180119T170449_0100.ZIP. I've attached an image with screenshots corresponding to step 5. Any pointers on where I going wrong are greatly appreciated. Thanks!
Mandy

Step 1
The tsm_nn netCDF that contains the TSM values but no latitude or longitude information. Per this product information https://earth.esa.int/eogateway/documents/20142/37627/MERIS-Sentinel-3-Like-L1-andL2-PFS.pdf I add the latitude and longitude from the geo_coordinates netCDF to the tsm_nn netCDF as spatial dimensions. This creates a new netCDF that has TSM along with latitude and longitude.

# NetCDF with data but missing lat/lon dimensions
tsm_nc = os.path.join(base_directory, "tsm_nn.nc")
# NetCDF containing lat/lon dimensions
geocoord_nc = os.path.join(base_directory, "geo_coordinates.nc")

if not os.path.exists(tsm_nc) or not os.path.exists(geocoord_nc):
raise FileNotFoundError("One or both NetCDF files are missing in the directory.")

ds_tsm = xr.open_dataset(tsm_nc)
ds_geocoord = xr.open_dataset(geocoord_nc)

# B) Extract latitude & longitude
lat = ds_geocoord["latitude"].values # Extract lat as NumPy array
lon = ds_geocoord["longitude"].values # Extract lon as NumPy array

# C) Ensure lat/lon are correctly formatted
if lat.ndim == 2: # If 2D, extract unique values along the correct axis
lat = lat[:, 0] # Take the first column (assuming lat is constant across rows)
if lon.ndim == 2:
lon = lon[0, :] # Take the first row (assuming lon is constant across columns)

print(f"Latitude shape after fix: {lat.shape}") # Should be (4289,)
print(f"Longitude shape after fix: {lon.shape}") # Should be (4481,)

# D) Extract Data & Ensure Correct Shape
var_name = list(ds_tsm.data_vars.keys())[0] # Get first variable name
data_values = ds_tsm[var_name].values # Extract data

# If data is 3D (e.g., time, lat, lon), select the first time step
if data_values.ndim == 3:
data_values = data_values[0, :, :]

# Ensure the data shape matches (lat, lon)
if data_values.shape != (len(lat), len(lon)):
raise ValueError(
f"Mismatch: data shape {data_values.shape} vs expected ({len(lat)}, {len(lon)})"
)

print(f"Final Data shape: {data_values.shape}") # Should match (4289, 4481)

# E) Create a new NetCDF dataset with correct dimensions
ds_new = xr.Dataset(
{
var_name: (["latitude", "longitude"], data_values)
},
coords={
"latitude": ("latitude", lat), # Assign dimensions explicitly
"longitude": ("longitude", lon)
}
)

# F) Save the updated dataset
tsm_nc_spdm = os.path.join(base_directory, "tsm_nn_spdm.nc")
ds_new.to_netcdf(tsm_nc_spdm)


Step 2
I convert this new netCDF into a geotiff. In the conversion I assign coordinate reference system EPSG 4326 to the geotiff based on this tutorial https://help.marine.copernicus.eu/en/articles/5029956-how-to-convert-netcdf-to-geotiff.

# Define path to input netCDF
tsm_nc_spdm = os.path.join(base_directory, "tsm_nn_spdm.nc")

# Open netcdf
tsm_nc_spdm = xr.open_dataset(tsm_nc_spdm)

# Extract variable
tsm = tsm_nc_spdm['TSM_NN']

# Define spatial dimentions
tsm = tsm.rio.set_spatial_dims(x_dim='longitude', y_dim='latitude')

# Assign CRS
tsm.rio.write_crs("epsg:4326", inplace=True)

# Define name of output GeoTIFF raster in the working directory
output_file = "tsm_nn_bad.tif"

# Export netcdf as geotiff raster
tsm.rio.to_raster(output_file)


Step 3
Attempt to plot the geotiff, which results in "ValueError: Axis limits cannot be NaN or Inf". Checking the geotiff data reveals bounding box and transformation are incorrect (Bounds: BoundingBox(left=nan, bottom=nan, right=nan, top=nan) Transform: | nan, nan, nan|| nan, nan, nan|| 0.00, 0.00, 1.00|)

def view_geotiff_with_basemap(file_path):
"""
Opens and displays a GeoTIFF file with a basemap, handling NaN/Inf values properly.

:param file_path: Path to the GeoTIFF file
"""
with rasterio.open(file_path) as src:
# Read data
data = src.read(1).astype(np.float32)

# Handle NoData values
if src.nodata is not None:
data[data == src.nodata] = np.nan

# Check if all values are NaN
if np.all(np.isnan(data)):
raise ValueError(f"Error: The dataset {file_path} contains only NaN values and cannot be plotted.")

# Compute extent
extent = [src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top]

# Mask invalid values
data = np.ma.masked_invalid(data)

# Create plot with Cartopy projection
fig, ax = plt.subplots(figsize=(10, 6), subplot_kw={"projection": ccrs.PlateCarree()})

# Add basemap features
ax.set_extent(extent, crs=ccrs.PlateCarree()) # Ensure correct extent
ax.add_feature(cfeature.LAND, edgecolor='black', facecolor='lightgray') # Land feature
ax.add_feature(cfeature.OCEAN, facecolor='lightblue') # Ocean feature
ax.add_feature(cfeature.BORDERS, linestyle=":") # Country borders
ax.add_feature(cfeature.COASTLINE, edgecolor="black") # Coastline
ax.gridlines(draw_labels=True, linestyle="--", color="gray")

# Plot the raster data
img = ax.imshow(data, cmap="viridis", extent=extent, origin="upper", transform=ccrs.PlateCarree())

# Add colorbar
cbar = plt.colorbar(img, ax=ax, orientation="vertical", fraction=0.03, pad=0.04)
cbar.set_label("Value")

# Set title
plt.title(f"GeoTIFF Visualization with Basemap: {file_path}")

# Show the plot
plt.show()

# Example Usage
tsm_tif_bad = "/Users/lopezama/Documents/Blackwood/MERIS/download_test/data/20120402T180502_use/tsm_nn_bad.tif"
view_geotiff_with_basemap(tsm_tif_bad)

# Check geotiff details to see why plot doesn't work
with rasterio.open("/Users/lopezama/Documents/Blackwood/MERIS/download_test/data/20120402T180502_use/tsm_nn_bad.tif") as src:
print("CRS:", src.crs) # CRS: None means missing spatial reference
print("Bounds:", src.bounds) # Bounds: [nan, nan, nan, nan] means corrupt transformation
print("Transform:", src.transform) # Transform: Affine(nan, nan, nan, nan, nan, nan means broken metadata
print("Data shape:", src.shape)


Step 4
Try to correct the bounding box and transformation. The ROI is the entire western US coastline from Washington State-Canada border to the California-Mexico border so the bounding box is large.

############### Step 4 #########################
# Fix geotiff bounding box and transformation

### ROI = entire western US coast so the bounding box is large

# Min lat (South): 32.4
# Max lat (North) 48.6
# Min lon (West): -125.2
# Max lon (East): -116.8

def fix_geotiff(file_path, output_path):
"""
Fixes the GeoTIFF file by manually setting the transform and CRS.

:param file_path: Path to the GeoTIFF file
:param output_path: Path to save the fixed GeoTIFF
"""
with rasterio.open(file_path) as src:
# Read the data
data = src.read(1).astype(np.float32)

# Define the bounding box (min_lon, max_lon, min_lat, max_lat)
min_lon = -125.2 # Set your correct min longitude
max_lon = -116.8 # Set your correct max longitude
min_lat = 32.4 # Set your correct min latitude
max_lat = 48.6 # Set your correct max latitude

# Define the transform (top-left corner, pixel size)
# from_origin(top_left_x, top_left_y, pixel_width, pixel_height)
transform = from_origin(min_lon, max_lat, (max_lon - min_lon) / data.shape[1], (max_lat - min_lat) / data.shape[0])

# Set the CRS (assuming WGS84 - EPSG:4326)
crs = "EPSG:4326"

# Write the fixed GeoTIFF
with rasterio.open(output_path, 'w', driver='GTiff',
count=1, dtype=data.dtype,
crs=crs, transform=transform,
width=data.shape[1], height=data.shape[0]) as dst:
dst.write(data, 1)

# Example Usage
tsm_tif_bad = "/Users/lopezama/Documents/Blackwood/MERIS/download_test/data/20120402T180502_use/tsm_nn_bad.tif"
tsm_tif_fix = "/Users/lopezama/Documents/Blackwood/MERIS/download_test/data/20120402T180502_use/tsm_nn_fix.tif"
fix_geotiff(tsm_tif_bad, tsm_tif_fix)

Step 5
Plot the modified geotiff. The bounding box properly shows the whole western US coastline, and the TSM data appear around southern California (which is correct for this specific image). But, some TSM values appear over land and far offshore in the Pacific Ocean. See attached image A = screenshot of whole MERIS swath, B = python plot from step 5, C = screenshot of geotiff from step 5 in QGIS with EPSG set to 4326.

# Plot modified geotiff with basemap

def view_geotiff_with_basemap(file_path):
"""
Opens and displays a GeoTIFF file with a basemap, handling NaN/Inf values properly.

:param file_path: Path to the GeoTIFF file
"""
with rasterio.open(file_path) as src:
# Read data
data = src.read(1).astype(np.float32)

# Handle NoData values
if src.nodata is not None:
data[data == src.nodata] = np.nan

# Check if all values are NaN
if np.all(np.isnan(data)):
raise ValueError(f"Error: The dataset {file_path} contains only NaN values and cannot be plotted.")

# Compute extent
extent = [src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top]

# Mask invalid values
data = np.ma.masked_invalid(data)

# Create plot with Cartopy projection
fig, ax = plt.subplots(figsize=(10, 6), subplot_kw={"projection": ccrs.PlateCarree()})

# Add basemap features
ax.set_extent(extent, crs=ccrs.PlateCarree()) # Ensure correct extent
ax.add_feature(cfeature.LAND, edgecolor='black', facecolor='lightgray') # Land feature
ax.add_feature(cfeature.OCEAN, facecolor='lightblue') # Ocean feature
ax.add_feature(cfeature.BORDERS, linestyle=":") # Country borders
ax.add_feature(cfeature.COASTLINE, edgecolor="black") # Coastline
ax.gridlines(draw_labels=True, linestyle="--", color="gray")

# Plot the raster data
img = ax.imshow(data, cmap="viridis", extent=extent, origin="upper", transform=ccrs.PlateCarree())

# Add colorbar
cbar = plt.colorbar(img, ax=ax, orientation="vertical", fraction=0.03, pad=0.04)
cbar.set_label("Value")

# Set title
plt.title(f"GeoTIFF Visualization with Basemap: {file_path}")

# Show the plot
plt.show()

# Example Usage
tsm_tif_fix = "/Users/lopezama/Documents/Blackwood/MERIS/download_test/data/20120402T180502_use/tsm_nn_fix.tif"
view_geotiff_with_basemap(tsm_tif_fix)

[OB.DAAC - SeanBailey]

I believe the issue may stem from the conversion to geotiff step. You are defining the CRS to be EPSG:4326. While the lat/lons in the L2 file are WGS84 coordinates, they are not in a regular gridded projection. This might be confusing the operator.

Sean

[lopezama]

Hi Sean, thanks for the response. I changed the approach based on your feedback, but I want to make sure I'm doing this properly. In this example I'm referring to the EN1_MDSI_MER_FRS_2P_2012040T184126_20120401T184409...ZIP shown in the attached screenshot. My code is below, but to summarize:

Steps 1-5 I assign the lat/lon WGS84 coordinates from the geo_coordinates.nc to the tsm_nn.nc as spatial dimensions, making a new TSM netCDF that has lat/lon info called tsm_nn_spdm.nc.

Step 6 I converted my ROI polygon shapefile from CRS EPSG 4326 to WGS84 coordinates.

Step 7 I plot the new TSM netCDF on a basemap with the ROI - see attached.

Question 1 - In the plot from step 7 there appear to be TSM values on land, which makes me think there is still something wrong with the georeferencing. Do you have suggestions on how to troubleshoot?

Question 2 - have the MERIS L2 TSM data have already gone through quality screening (i.e., cloud masks, land masks, etc. have been applied and errant data removed)? My understanding is that the L2 data are already screened based on 2.7.1.2.6.4 MERIS Ocean Colour Processing (step 2.9) https://earth.esa.int/eogateway/documents/20142/37627/MERIS-product-handbook.pdf, but I may be misunderstanding.

Thanks, Mandy

# Step 1: Load NetCDFs
tsm_nc = os.path.join(base_directory, "tsm_nn.nc")
geocoord_nc = os.path.join(base_directory, "geo_coordinates.nc")

if not os.path.exists(tsm_nc) or not os.path.exists(geocoord_nc):
raise FileNotFoundError("One or both NetCDF files are missing in the directory.")

ds_tsm = xr.open_dataset(tsm_nc)
ds_geocoord = xr.open_dataset(geocoord_nc)

# Step 2: Extract latitude & longitude
lat = ds_geocoord["latitude"].values
lon = ds_geocoord["longitude"].values

# Step 3: Extract TSM data and ensure correct shape
var_name = list(ds_tsm.data_vars.keys())[0] # Get the first variable name
data_values = ds_tsm[var_name].values

# Step 3A: If data is 3D (e.g., time, lat, lon), select the first time step
if data_values.ndim == 3:
data_values = data_values[0, :, :]

# Step 3B: Expected shape from TSM data
expected_shape = data_values.shape

# Step 3C: Debugging: Print original shapes
print(f"Original Shapes:")
print(f" - Latitude: {lat.shape}")
print(f" - Longitude: {lon.shape}")
print(f" - Data ({var_name}): {data_values.shape}")

# Step 3D: If lat/lon are 1D, convert to 2D
if lat.ndim == 1 and lon.ndim == 1:
lon, lat = np.meshgrid(lon, lat)

# Step 3E: Debugging: Print shapes after conversion
print(f"Shapes After Conversion (if applied):")
print(f" - Latitude: {lat.shape}")
print(f" - Longitude: {lon.shape}")
print(f" - Data ({var_name}): {data_values.shape}")

# Step 3F: Ensure final shapes match
if lat.shape != expected_shape or lon.shape != expected_shape:
raise ValueError(f"Lat/Lon shape mismatch: Expected {expected_shape}, Got {lat.shape} & {lon.shape}")

# Step 4: Create new NetCDF dataset
ds_new = xr.Dataset(
{
var_name: (["y", "x"], data_values)
},
coords={
"latitude": (["y", "x"], lat),
"longitude": (["y", "x"], lon)
}
)

# Step 5: Save NetCDF
output_nc = os.path.join(base_directory, "tsm_nn_spdm_test.nc")
ds_new.to_netcdf(output_nc)

# Step 6
def convert_shapefile_to_latlon(input_shapefile, output_shapefile):
"""
Converts a shapefile's CRS to latitude and longitude (EPSG:4326).

:param input_shapefile: Path to the input shapefile.
:param output_shapefile: Path to save the converted shapefile.
"""
# Load the shapefile
gdf = gpd.read_file(input_shapefile)

# Check current CRS
if gdf.crs is None:
raise ValueError("Error: Input shapefile has no CRS. Please define it first.")

# Convert to latitude and longitude (WGS 84)
gdf = gdf.to_crs(epsg=4326)

# Save the transformed shapefile
gdf.to_file(output_shapefile, driver="ESRI Shapefile")

print(f"Shapefile successfully converted to latitude/longitude and saved to: {output_shapefile}")