Introduction and Background
Geographical information systems (GIS) often involve working with spatial data, including raster files that represent data on a grid. NetCDF (Network Common Data Form) is a popular format for storing multidimensional scientific data, including raster datasets. In this tutorial, we’ll explore how to extract geographical coordinates from a NetCDF raster file using Python.
Following dataset is used in the tutorial:
NetCDF dataset
1. NetCDF Format
NetCDF is a self-describing, machine-independent data format for representing scientific data. It is commonly used for climate and weather data, satellite imagery, and various other geospatial datasets. NetCDF files can store multidimensional data, making them well-suited for applications where data has dimensions such as time, latitude, and longitude.
2. Python Libraries Used
- xarray: Xarray is a powerful library designed for working with labeled multidimensional arrays, particularly NetCDF files. It simplifies data manipulation and analysis, providing a high-level interface for working with labeled data.
- numpy: Numpy is a fundamental package for scientific computing in Python. It provides support for large, multi-dimensional arrays and matrices, along with mathematical functions to operate on these arrays.
- pandas: Pandas is a data manipulation and analysis library. It provides data structures such as DataFrame, which is ideal for organizing and analyzing structured data.
- geopandas: Geopandas is geospatial manipulation library for vector data processing. Here we are using it to convert raster coordinates into shapefile
- rioxarray: Rioxarray extends xarray’s capabilities by providing geospatial functionality, such as reprojecting, getting coordinate reference system (CRS) information, and handling spatial data.
3. Tutorial Steps
Step 1: Install Required Libraries
|
1 2 3 |
(base) geoknight@pop-os:~$conda create -n spatial-dev.guru python=3.11 (base) geoknight@pop-os:~$conda activate spatial-dev.guru (spatial-dev.guru) geoknight@pop-os:~$conda install -c conda-forge xarray matplotlib numpy pandas rioxarray netcdf geopandas |
This command installs the necessary Python libraries (xarray, matplotlib, numpy, pandas, and rioxarray) using the conda
Step 2: Import Libraries
|
1 2 3 4 5 |
import xarray as xr import numpy as np import pandas as pd import geopandas as gpd import math |
These lines import the required libraries for working with NetCDF files (xarray), numerical operations (numpy), and data manipulation (pandas and geopandas).
Step 3: Read NetCDF File Using xarray
ds = xr.open_dataset("HLSTimeSeries.nc", decode_coords="all")Here, xr.open_dataset is used to open the NetCDF file named “HLSTimeSeries.nc”. The decode_coords="all" argument ensures that coordinate decoding is performed for compatibility with rioxarray.
Step 4: Check Coordinate Reference System (CRS) of NetCDF
ds.rio.crsThis line checks and prints the Coordinate Reference System (CRS) of the NetCDF file using the rioxarray extension of xarray.
Step 5: Reproject NetCDF to EPSG:4326
ds_4326 = ds.rio.reproject("EPSG:4326")This code snippet reprojects the NetCDF file to EPSG:4326, which is a common geographic coordinate reference system.
Step 6: Get XY Coordinates in 32756 Projection
x, y = np.meshgrid(ds.x, ds.y)This section creates a meshgrid of XY coordinates in the original projection (32756) using NumPy’s meshgrid function.
Step 7: Get All Variables in the NetCDF File
variables = list(ds.var())These lines list all the variables present in the NetCDF file.
Step 8: Convert Variables and Bands to DataFrame
|
1 2 3 4 5 6 7 8 |
variable_name = variables[0] raster_values = ds[variable_name][0].values df = pd.DataFrame({ "x": x.flatten(), "y": y.flatten(), "raster_value": raster_values.flatten() }) |
This section selects the first variable, extracts the raster values for the first band, and creates a pandas DataFrame (df) containing the x, y coordinates, and raster values.
Step 9: Convert Dataframe to GeoDataFrame
|
1 2 |
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['x'], df['y']), crs=ds.rio.crs) gdf[~gdf['raster_value'].apply(lambda x: math.isnan(x))].to_file("shp/points.shp") |
Feel free to customize the code and explanations based on your specific needs and dataset.
Complete Code:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 |
# Import Libraries import xarray as xr import numpy as np import pandas as pd import geopandas as gpd import math # Use xarray to read netcdf file. It is recommended to use decode_coords="all" for compatibility with RioXarray. You can rioxarray functionalities like reprojecting, getting CRS, extent etc information ds = xr.open_dataset("HLSTimeSeries.nc", decode_coords="all") # Check CRS of netcdf ds.rio.crs float(ds.x[0]), float(ds.y[0]) # Reprojecting netcdf to 4326 ds_4326 = ds.rio.reproject("EPSG:4326") float(ds_4326.x[0]), float(ds_4326.y[0]) # Get XY coordinates which are in 32756 projection. If you want Geographic coordnates, then you can reproject the xarray ataset to 4326 x, y = np.meshgrid(ds.x, ds.y) # Get all variables variables = list(ds.var()) # Converting a variables and its band into dataframe of xy coordinates and its raster values raster_values = ds[variables[0]][0].values df = pd.DataFrame({"x": x.flatten(), "y": y.flatten(), "raster_value": raster_values.flatten()}) # Converting to GeoDataFrame gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy( df['x'], df['y']), crs=ds.rio.crs) gdf[~gdf['raster_value'].apply( lambda x: math.isnan(x))].to_file("shp/points.shp") |
I hope this tutorial will create a good foundation for you. If you want tutorials on another GIS topic or you have any queries, please send an email at contact@spatial-dev.guru.