Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixing the GOES transformation to regular lat/lon COG #1

Open
j08lue opened this issue Oct 24, 2024 · 0 comments
Open

Fixing the GOES transformation to regular lat/lon COG #1

j08lue opened this issue Oct 24, 2024 · 0 comments

Comments

@j08lue
Copy link

j08lue commented Oct 24, 2024

So the conversion of GOES from NetCDF to COG did not produce a correctly geo-referenced GeoTIFF with

https://github.com/NASA-IMPACT/veda-cyclone/blob/56037d536aa801fa3f6a767b28e9c8c53e80b4de/transformation_script/convert_netcdf_cog_GOES.ipynb

The issue was that rioxarray/rasterio/GDAL does not work well with 2D coordinate variables created from transforming the GOES coordinates to lat/lon.

Instead, you need to reproject the data to a regular lat/lon grid.

This works and created correctly geo-referenced COGs:

from pyproj import CRS

def goes_to_wgs(ds, variable_name):
    
    # convert coordinates to meters https://gis.stackexchange.com/a/350006
    sat_height = ds["goes_imager_projection"].attrs["perspective_point_height"]
    ds = ds.assign_coords({
        "x": ds["x"].values * sat_height,
        "y": ds["y"].values * sat_height,
    })
    
    crs = CRS.from_cf(ds["goes_imager_projection"].attrs)
    ds.rio.write_crs(crs.to_string(), inplace=True)
    
    da = ds[variable_name]

    return da.rio.reproject("epsg:4326")

Create a COG from a NetCDF file on disk:

curl -O https://noaa-goes16.s3.amazonaws.com/ABI-L1b-RadF/2017/129/11/OR_ABI-L1b-RadF-M3C03_G16_s20171291115392_e20171291126159_c20171291126196.nc
import xarray as xr
import rioxarray

path = "OR_ABI-L1b-RadF-M3C03_G16_s20171291115392_e20171291126159_c20171291126196.nc"

with xr.open_dataset(path, engine='h5netcdf') as ds:
    da = goes_to_wgs(ds, variable_name="Rad")
    COG_PROFILE = {"driver": "COG", "compress": "DEFLATE", "predictor": 2}
    da.rio.to_raster("goes16.tif", **COG_PROFILE)

or straight from AWS S3:

import s3fs 
import xarray as xr
import rioxarray

fs = s3fs.S3FileSystem(anon=True)

path = "s3://noaa-goes16/ABI-L1b-RadF/2017/129/11/OR_ABI-L1b-RadF-M3C03_G16_s20171291115392_e20171291126159_c20171291126196.nc"

with fs.open(path) as fileObj:
    with xr.open_dataset(fileObj, engine='h5netcdf') as ds:
        da = goes_to_wgs(ds, variable_name="Rad")
        COG_PROFILE = {"driver": "COG", "compress": "DEFLATE", "predictor": 2}
        da.rio.to_raster("goes16.tif", **COG_PROFILE)

Notebook here: https://gist.github.com/j08lue/f949404cf95c9a521b6b72f33a4cbf3b

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant