Pyproj wraps the proj4 C library, which is very widely used for coordinate system transformations. If you have an XY pair in projected coordinates and you know the definition of the projected coordinate system, you can use pyproj's transform function like this:
import pyproj
# Output coordinates are in WGS 84 longitude and latitude
projOut = pyproj.Proj(init = 'epsg:4326')
# Input coordinates are in meters on the Polar Stereographic
# projection given in the netCDF file
projIn = pyproj.Proj('+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +
k = 1 + x_0 = 0 + y_0 = 0 + datum = WGS84 + units = m + no_defs ',
preserve_units = True)
# here is a coordinate pair near the middle of your data set
x, y = 0.0, -2000000
# transform x, y to lon / lat
lon, lat = pyproj.transform(projIn, projOut, x, y)
# answer: lon = -45.0;
lat = 71.6886Setting projIn is harder, because your netCDF file defines its coordinate system with a WKT string, which (I'm pretty sure) can't be read directly by proj4 or pyproj. However, pyproj.Proj() will take a proj4 parameter string as an argument. I've already given you the one you need for this operation, so you can just take my for for it that this
+proj = stere + lat_0 = 90 + lat_ts = 70 + lon_0 = -45 + k = 1 + x_0 = 0 + y_0 = 0 + datum = WGS84 + units = m + no_defs
is the equivalent of this (which is copied directly from your netCDF file):
PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84", 6378137, 298.257223563,
AUTHORITY["EPSG", "7030"]],
AUTHORITY["EPSG", "6326"]],
PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]],
UNIT["degree", 0.0174532925199433, AUTHORITY["EPSG", "9122"]],
AUTHORITY["EPSG", "4326"]],
PROJECTION["Polar_Stereographic"],
PARAMETER["latitude_of_origin", 70],
PARAMETER["central_meridian", -45],
PARAMETER["scale_factor", 1],
PARAMETER["false_easting", 0],
PARAMETER["false_northing", 0],
UNIT["metre", 1, AUTHORITY["EPSG", "9001"]],
AXIS["X", EAST],
AXIS["Y", NORTH],
AUTHORITY["EPSG", "3413"]]
'FYI, you can use file.nc?varname to access a single NetCDF variable as per https://docs.generic-mapping-tools.org/6.2/cookbook/features.html#modifiers-for-coards-compliant-netcdf-files. E.g. gmt grdinfo BedMachineAntarctica_2020-07-15_v02.nc?bed will give:,I tried running gmt grdproject on BedMachine Antarctica but it took a long time and didn’t finish. If you have a smaller study region (e.g. a single glacier/ice stream), it might work though. Or maybe just stick with mapproject ,Thanks for the tip on how to select one variable from a netCDF file. Looks like I’m not very good at digging into the documentation (!) Regards, David,And just to busybody a bit, is there a reason to convert the BedMachine Antarctica NetCDF grid to an XYZ file before reprojecting? Just wondering if using grdproject might be better, though I’m not sure what you’re trying to do exactly.
Try using the EPSG code (3031 for Antarctic Polar Stereographic). The below seems to give sensible-ish results, I’m not sure why the mapproject -E option doesn’t work as expected with the inverse -I option.
gmt mapproject XYZ.txt - JEPSG: 3031 - R - 180 / 180 / -90 / -60 - I > LONLATZ.txt
produces the following (note that it’s Longitude, Latitude, Elevation). Not sure if the Elevation is correct, but at least it looks better than before…
-45 - 48.4643831064 - 5915.61719 - 45.0430083837 - 48.4929562939 - 5962.10742 - 45.0860813347 - 48.5215132292 - 5959.86377 - 45.1292189499 - 48.5500538498 - 5866.03711 - 45.1724213258 - 48.5785780935 - 5661.64453 - 45.2156885592 - 48.6070858975 - 5622.0415 - 45.2590207468 - 48.6355771991 - 5500.5874 - 45.3024179853 - 48.6640519356 - 5393.52588 - 45.3458803713 - 48.6925100438 - 5370.03662 - 45.3894080017 - 48.7209514607 - 5484.46045 - 45.4330009729 - 48.7493761231 - 5668.95605 - 45.4766593818 - 48.7777839675 - 5865.19629 - 45.5203833248 - 48.8061749305 - 5878.93457 - 45.5641728987 - 48.8345489484 - 5886.30225 - 45.6080282001 - 48.8629059574 - 5934.98633 - 45.6519493255 - 48.8912458938 - 5999.31201 - 45.6959363715 - 48.9195686934 - 6071.00391 - 45.7399894347 - 48.9478742921 - 6139.77783 - 45.7841086116 - 48.9761626256 - 6057.91699 - 45.8282939987 - 49.0044336294 - 6030.22168
FYI, you can use file.nc?varname to access a single NetCDF variable as per https://docs.generic-mapping-tools.org/6.2/cookbook/features.html#modifiers-for-coards-compliant-netcdf-files. E.g. gmt grdinfo BedMachineAntarctica_2020-07-15_v02.nc?bed will give:
BedMachineAntarctica_2020 - 07 - 15_ v02.nc: Title: bed topography
BedMachineAntarctica_2020 - 07 - 15_ v02.nc: Command:
BedMachineAntarctica_2020 - 07 - 15_ v02.nc: Remark:
BedMachineAntarctica_2020 - 07 - 15_ v02.nc: Gridline node registration used[Cartesian grid]
BedMachineAntarctica_2020 - 07 - 15_ v02.nc: Grid file format: nf = GMT netCDF format(32 - bit float), CF - 1.7
BedMachineAntarctica_2020 - 07 - 15_ v02.nc: x_min: -3333000 x_max: 3333000 x_inc: 500 name: Cartesian x - coordinate[meter] n_columns: 13333
BedMachineAntarctica_2020 - 07 - 15_ v02.nc: y_min: -3333000 y_max: 3333000 y_inc: 500 name: Cartesian y - coordinate[meter] n_rows: 13333
BedMachineAntarctica_2020 - 07 - 15_ v02.nc: v_min: 0 v_max: 0 name: bed topography[meters]
BedMachineAntarctica_2020 - 07 - 15_ v02.nc: scale_factor: 1 add_offset: 0
BedMachineAntarctica_2020 - 07 - 15_ v02.nc: format: netCDF - 4 chunk_size: 953, 953 shuffle: on deflation_level: 1The netCDF file is not following the CF-convention, I would say it follows none. CDO is following the CF-convention and needs dimensions like time, lat, lon with the appropriate units. And if the dimension names are not lat and lon the variables have to have a coordinates attribute which leads CDO to the grid dimensions. Furthermore, the variable dimensions should be (time,lat,lon) or (time,lev,lat,lon).,I have try to use CDO functions like -setattribute or remapbil, trying to assign coordinate units to Lon/Lat variables or to interpolate from another NetCDF with the similar grid, but I have no results.,Unfortunately, the original NetCDF file does not follow the conventions you indicate, they are outputs of a hydrological model that I do not have access to, so I must modify the file on my own, in order to process and use the information in the file.,You have to be careful because the variable ET (on horizontal grid) has just one time step which refers to the original time variable dim_time_2, which is missing. The variables agno and mes (both are time series) refer to the original dim_time variable, which is also missing. CDO accepts only one time dimension that means that you should seperate the variables in two different files.
This is the information of my NetCDF file:
File format: NetCDF - 1: Institut Source T Steptype Levels Num Points Num Dtype: Parameter name 1: unknown unknown c instant 1 1 400 1 F64: lat 2: unknown unknown c instant 1 1 220 2 F64: lon 3: unknown unknown c instant 1 1 444 3 F64: agno 4: unknown unknown c instant 1 1 444 3 F64: mes 5: unknown unknown c instant 444 2 88000 4 F32: ET 6: unknown unknown c instant 444 2 88000 4 F32: runoff 7: unknown unknown c instant 444 2 88000 4 F32: fsca 8: unknown unknown c instant 444 2 88000 4 F32: riego 9: unknown unknown c instant 444 2 88000 4 F32: SWE 10: unknown unknown c instant 444 2 88000 4 F32: SM 11: unknown unknown c instant 444 2 88000 4 F32: recarga Grid coordinates: 1: generic: points = 400 2: generic: points = 220 3: generic: points = 444 4: generic: points = 88000(400 x220) Vertical coordinates: 1: surface: levels = 1 2: generic: levels = 444
I have already use the "setattribute" operator, and a I get something like this:
File format: NetCDF - 1: Institut Source T Steptype Levels Num Points Num Dtype: Parameter name 1: unknown unknown c instant 1 1 444 1 F64: agno 2: unknown unknown c instant 1 1 444 1 F64: mes 3: unknown unknown c instant 1 2 88000 2 F32: ET Grid coordinates: 1: generic: points = 444 2: generic: points = 88000(400 x220) Vertical coordinates: 1: surface: levels = 1 2: generic: levels = 1 cdo sinfon: Processed 3 variables(0.03 s 57 MB)ncdump -h regionalizacion_1979_2015_ET_sel.nc
netcdf regionalizacion_1979_2015_ET_sel {
dimensions: dim_lat = 400;
dim_lon = 220;
dim_time = 444;
dim_time_2 = 1;
variables: double lat(dim_lat);
lat: Latitud = "norte-sur";
double lon(dim_lon);
lon: Longitud = "oeste-este";
double agno(dim_time);
agno: inicio = "1979/01";
agno: fin = "2015/12";
double mes(dim_time);
float ET(dim_time_2, dim_lon, dim_lat);
ET: Long_name = "ET_sin_riego";
ET: unidades = "[mm/d]";
// global attributes:
: CDI = "Climate Data Interface version 1.9.2 (http://mpimet.mpg.de/cdi)";: Conventions = "CF-1.6";: history = "Mon Jan 06 15:41:25 2020: cdo select,name=lat,lon,agno,mes,ET,level=1 regionalizacion_1979_2015.nc regionalizacion_1979_2015_ET_sel.nc";: CDO = "Climate Data Operators version 1.9.2 (http://mpimet.mpg.de/cdo)";
}1. Reorder the dimensions of ET with NCO's ncpdq
ncpdq - O - a dim_time_2, dim_lat, dim_lon $infile tmp.nc
Result:
dimensions:
dim_lat = 400;
dim_lon = 220;
dim_time = 444;
dim_time_2 = 1;
variables:
double lat(dim_lat);
lat: Latitud = "norte-sur";
double lon(dim_lon);
lon: Longitud = "oeste-este";
double agno(dim_time);
agno: inicio = "1979/01";
agno: fin = "2015/12";
double mes(dim_time);
float ET(dim_time_2, dim_lat, dim_lon);
ET: Long_name = "ET_sin_riego";
ET: unidades = "[mm/d]";3. Change the set of attributes for ET, lat, and lon.
cdo - setattribute, FILE = attr_file.txt tmp.nc tmp2.nc
4. Select the correct grid, sellonlatbox, and change the time dimension name. In the original netCDF file
the dim_time_2 variable doesn't exist, that's why it's deleted, here.
cdo--reduce_dim - sellonlatbox, -75, -70, -30, -20 - selgrid, 2 tmp2.nc tmp3.nc26 August 2013
Here's some output, relevant to the problem we'll present, from "ncdump -h" applied to the example file:
dimensions:
MT = UNLIMITED; // (1 currently) Y = 850 ; X = 712 ; Depth = 10 ;variables: double MT(MT) ; MT:units = "days since 1900-12-31 00:00:00" ; double Date(MT) ; Date:units = "day as %Y%m%d.%f" ; float Depth(Depth) ; Depth:units = "m" ; Depth:positive = "down" ; int Y(Y) ; Y:axis = "Y" ; int X(X) ; X:axis = "X" ; float Latitude(Y, X) ; Latitude:units = "degrees_north" ; float Longitude(Y, X) ; Longitude:units = "degrees_east" ; float temperature(MT, Depth, Y, X) ; temperature:coordinates = "Longitude Latitude Date" ; temperature:standard_name = "sea_water_potential_temperature" ; temperature:units = "degC" ; float salinity(MT, Depth, Y, X) ; salinity:coordinates = "Longitude Latitude Date" ; salinity:standard_name = "sea_water_salinity" ; salinity:units = "psu" ; MT = UNLIMITED; // (1 currently) Y = 850 ; X = 712 ; Depth = 10 ;variables: double MT(MT) ; MT:units = "days since 1900-12-31 00:00:00" ; double Date(MT) ; Date:units = "day as %Y%m%d.%f" ; float Depth(Depth) ; Depth:units = "m" ; Depth:positive = "down" ; int Y(Y) ; Y:axis = "Y" ; int X(X) ; X:axis = "X" ; float Latitude(Y, X) ; Latitude:units = "degrees_north" ; float Longitude(Y, X) ; Longitude:units = "degrees_east" ; float temperature(MT, Depth, Y, X) ; temperature:coordinates = "Longitude Latitude Date" ; temperature:standard_name = "sea_water_potential_temperature" ; temperature:units = "degC" ; float salinity(MT, Depth, Y, X) ; salinity:coordinates = "Longitude Latitude Date" ; salinity:standard_name = "sea_water_salinity" ; salinity:units = "psu" ;A naive metric for distance squared between (lat,lon) and (lat0, lon0) just uses the "Euclidean norm", because it's easy to compute:
dist_squared = (lat - lat0) 2 + (lon - lon0) 2
MT = UNLIMITED; // (1 currently) Y = 850 ; X = 712 ; Depth = 10 ;variables: double MT(MT) ; MT:units = "days since 1900-12-31 00:00:00" ; double Date(MT) ; Date:units = "day as %Y%m%d.%f" ; float Depth(Depth) ; Depth:units = "m" ; Depth:positive = "down" ; int Y(Y) ; Y:axis = "Y" ; int X(X) ; X:axis = "X" ; float Latitude(Y, X) ; Latitude:units = "degrees_north" ; float Longitude(Y, X) ; Longitude:units = "degrees_east" ; float temperature(MT, Depth, Y, X) ; temperature:coordinates = "Longitude Latitude Date" ; temperature:standard_name = "sea_water_potential_temperature" ; temperature:units = "degC" ; float salinity(MT, Depth, Y, X) ; salinity:coordinates = "Longitude Latitude Date" ; salinity:standard_name = "sea_water_salinity" ; salinity:units = "psu" ;tile is an array of integers from 0 to 12, one for each tile of the lat-lon-cap grid.,k is an array of integers from 0 to 49 indicating the z_grid_index along this tile’s Z axis.,j is an array of integers from 0 to 89 indicating the y_grid_index along this tile’s Y axis.,i is an array of integers from 0 to 89 indicating the x_grid_index along this tile’s X axis.
[1]:import numpy as np
import xarray as xr
import sys
import matplotlib.pyplot as plt %
matplotlib inline
import json[2]:# # Import the ecco_v4_py library into Python
# # === === === === === === === === === === === === === ==
# #--If ecco_v4_py is not installed in your local Python library,
# # tell Python where to find it.For example,
if your ecco_v4_py
# # files are in /Users/ifenty / ECCOv4 - py / ecco_v4_py, then use:
sys.path.append('/home/ifenty/ECCOv4-py')
import ecco_v4_py as ecco[3]:# # Set top - level file directory for the ECCO NetCDF files # # === === === === === === === === === === === === === === === === === === === === === == # base_dir = '/home/username/' base_dir = '/home/ifenty/ECCOv4-release' # # define a high - level directory for ECCO fields ECCO_dir = base_dir + '/Release3_alt'