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)