geoscience-data

# Geoscience data Gridded and geospatial data carry their meaning in metadata: units, fill values, calendars, the coordinate reference system and the cell geometry. Most wrong answers come from ignoring one of them, such as averaging over a latitude-longitude grid without weights, treating `_FillValue` as zero, or mixing degrees with metres. ## NetCDF and HDF5 1. `xr.open_dataset(path)` decodes CF metadata: `scale_factor` and `add_offset` unpack integers, `_FillValue` and `missing_value` become NaN, and time is decoded from `units` ("days since 1850-01-01") with the `calendar` attribute. Use `decode_times=False` to see raw numbers, `use_cftime=True` for `noleap` or `360_day` calendars and for dates outside the range of nanosecond `datetime64`. Multi-file: `xr.open_mfdataset(pattern, combine="by_coords", chunks={...})` with dask. Inspect with `ncdump -h` or `ds.info()`, and read `units`, `standard_name`, `cell_methods` and `bounds` for every variable you use. HDF5 that is not NetCDF-4 opens with `h5py`. 2. Dimensions are axes (`time`, `lat`, `lon`); coordinates are labels on them. On a curvilinear grid `lat(y, x)` is a two-dimensional coordinate and `sel` by latitude no longer works directly. Longitude runs 0 to 360 in many models and -180 to 180 in others; convert with `ds.assign_coords(lon=((ds.lon + 180) % 360) - 180).sortby("lon")`. Latitude often descends from 90 to -90, so `sel(lat=slice(a, b))` needs that order. 3. Time stamps of monthly means may mark the start, middle or end of the month; consult `time_bnds`. Yearly means from monthly means must weight by `ds.time.dt.days_in_month`. `ds.resample(time="1MS").mean()` for calendar months. Times are UTC unless the file says otherwise. 4. Writing: `ds.to_netcdf(path, encoding=enc)` with a per-variable `enc` such as `{"var": {"dtype": "float32", "_FillValue": -9999.0, "zlib": True, "complevel": 4}}`; keep `units` and `standard_name` attributes on every variable, and never let a fill value coincide with a valid data value. ## Area weighting and regridding 5. On a regular latitude-longitude grid cell area scales with cos(latitude): `w = np.cos(np.deg2rad(ds.lat)); ds.weighted(w).mean(("lat", "lon"))`. Prefer an explicit `cell_area` variable or areas computed from `lat_bnds` when present. Unweighted means over-represent the poles. Totals (mass, volume, emissions) need true areas in square metres and a stated Earth radius (6371 km) or ellipsoid. 6. Regridding: conservative for fluxes and extensive quantities, bilinear for smooth intensive fields, nearest for categorical data (xesmf implements all three; `xr.interp` is bilinear). Name the method and the target grid. ## Raster 7. `rasterio.open(path)` exposes `crs`, `transform` (affine pixel-to-world), `nodata`, `res`, `bounds` and `tags()`. `src.read(1, masked=True)` honours nodata. GeoTIFF `AREA_OR_POINT` says whether a pixel value describes the cell or its centre. Reproject with `rasterio.warp.calculate_default_transform` and `reproject(..., resampling=Resampling.bilinear)`; `Resampling.nearest` for classes, `average` or `sum` when coarsening. GDAL equivalents: `gdalinfo -stats`, `gdalwarp -t_srs EPSG:3035 -tr 100 100 -r bilinear -dstnodata -9999 in.tif out.tif`, `gdal_translate -of COG`. `rioxarray.open_rasterio(path)` gives an xarray with `.rio.crs`, `.rio.reproject("EPSG:...")` and `.rio.write_nodata(...)`. ## Coordinate reference systems 8. EPSG:4326 is WGS 84 latitude-longitude in degrees; its authority axis order is latitude first while most software expects longitude first, so create transformers with `pyproj.Transformer.from_crs(4326, 3857, always_xy=True)`. EPSG:3857 (Web Mercator) distorts area and is for tiles, not measurement. UTM zones are EPSG:326xx north and 327xx south; equal-area choices include EPSG:3035 (Europe LAEA), EPSG:5070 (CONUS Albers) and EPSG:6933 (global EASE-Grid 2.0). `pyproj.CRS.from_epsg(code)` exposes `axis_info`, `is_geographic` and units. 9. Distances and areas on geographic coordinates are computed geodesically, not in degrees: `pyproj.Geod(ellps="WGS84").inv(lon1, lat1, lon2, lat2)` returns metres, and `geod.geometry_area_perimeter(polygon)` returns square metres. Haversine is spherical and off by up to half a percent; say which formula was used. Datum shifts (WGS 84, NAD83, ETRS89) matter at the metre level; elevation needs a stated vertical datum. ## Vector 10. `geopandas.read_file(path)` for shapefile, GeoPackage and GeoJSON; `gdf.crs`, `gdf.to_crs(epsg=...)`, `gdf.estimate_utm_crs()` for a local metric CRS, then `.area` and `.length` in that CRS's units (geographic CRS results are warned about and wrong). `sjoin(predicate="intersects")`, `overlay`, `dissolve`; validate with `is_valid` and `make_valid`. Shapefiles limit field names to 10 characters, files to 2 GB, one geometry type per layer and need `.shp`, `.shx`, `.dbf` and `.prj` together; prefer GeoPackage. Buffers in degrees are meaningless; buffer in a projected CRS. ## Units and reporting - Convert explicitly (K to degrees Celsius, kg m^-2 s^-1 to mm/day by 86400, Pa to hPa) or with `pint` and `cf_xarray`; put units in every column header. - State CRS (EPSG), resolution, resampling, weighting, calendar, fill handling, and the provenance of each file (URL, version, download date). ## Sources - CF Metadata Conventions: https://cfconventions.org/ - xarray user guide (I/O, weighted reductions, time series): https://docs.xarray.dev/en/stable/user-guide/ - netCDF4-python: https://unidata.github.io/netcdf4-python/ - rasterio reprojection: https://rasterio.readthedocs.io/en/stable/topics/reproject.html - GDAL `gdalwarp`: https://gdal.org/en/stable/programs/gdalwarp.html - pyproj `Transformer` and `Geod`: https://pyproj4.github.io/pyproj/stable/api/ - GeoPandas projections guide: https://geopandas.org/en/stable/docs/user_guide/projections.html - EPSG Geodetic Parameter Registry: https://epsg.org/