Subset NETCDF file with similar gdalinfo

[danhamill]

I am trying to create a georeferrenced NETCDF file using a NCAR temperature dataset available from the NCAR Data Gateway. The data files are served in lat/lon format and are not georeferenced. From the paper, they indicate the domain is based upon a DEM elevation dataset resampled to a 2km resolution. In my scripts I call it pr_targ.

Here is what the source data looks as an xarrary dataset:

<xarray.Dataset>
Dimensions:    (lat: 925, lon: 1280, time: 365)
Coordinates:
  * time       (time) datetime64[ns] 1981-01-01 1981-01-02 ... 1981-12-31
Dimensions without coordinates: lat, lon
Data variables:
    corr_pcp   (time, lat, lon) float32 ...
    elevation  (lat, lon) float64 ...
    latitude   (lat, lon) float64 ...
    longitude  (lat, lon) float64 ...
    pcp        (time, lat, lon) float32 ...
    snow_frac  (time, lat, lon) float64 ...
    t_mean     (time, lat, lon) float32 ...
    t_range    (time, lat, lon) float32 ...
    wetPert    (lat, lon) float64 ...
    windPert   (lat, lon) float64 ...
Attributes:
    NCO:          netCDF Operators version 4.7.9 (Homepage = http://nco.sf.ne...
    history:      Version 1.0 of the NCAR Alaska ensemble meteorological data...
    institution:  National Center for Atmospheric Research (NCAR), Boulder, C...
    title:        Ensemble (2km) Gridded Daily Temperature and Bias Adjusted ...
    source:       Version 1.0 of the time hierarchical GMET system using TIER...
    references:   Newman et al. 2020: Probabilistic Spatial Meteorological Es...

Add corresponding gdalinfo

C:\Users\Public\Desktop\OSGeo4W>gdalinfo "G:\Newman_AK_Data\alaska_daily_1981_050.nc4"
Warning 1: Recode from UTF-8 to CP_ACP failed with the error: "Invalid argument".
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Driver: netCDF/Network Common Data Format
Files: G:\Newman_AK_Data\alaska_daily_1981_050.nc4
Size is 512, 512
Metadata:
  NC_GLOBAL#history=Version 1.0 of the NCAR Alaska ensemble meteorological dataset, created November 2019
  NC_GLOBAL#institution=National Center for Atmospheric Research (NCAR), Boulder, CO USA
  NC_GLOBAL#NCO=netCDF Operators version 4.7.9 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)
  NC_GLOBAL#references=Newman et al. 2020: Probabilistic Spatial Meteorological Estimates for Alaska and the Yukon, JGR-Atmospheres
  NC_GLOBAL#source=Version 1.0 of the time hierarchical GMET system using TIER v1.0 climatological fields
  NC_GLOBAL#title=Ensemble (2km) Gridded Daily Temperature and Bias Adjusted Precipitation for Alaska and the Yukon (1980-2013)
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"G:\Newman_AK_Data\alaska_daily_1981_050.nc4":corr_pcp
  SUBDATASET_1_DESC=[365x925x1280] corr_pcp (32-bit floating-point)
  SUBDATASET_2_NAME=NETCDF:"G:\Newman_AK_Data\alaska_daily_1981_050.nc4":elevation
  SUBDATASET_2_DESC=[925x1280] elevation (64-bit floating-point)
  SUBDATASET_3_NAME=NETCDF:"G:\Newman_AK_Data\alaska_daily_1981_050.nc4":latitude
  SUBDATASET_3_DESC=[925x1280] latitude (64-bit floating-point)
  SUBDATASET_4_NAME=NETCDF:"G:\Newman_AK_Data\alaska_daily_1981_050.nc4":longitude
  SUBDATASET_4_DESC=[925x1280] longitude (64-bit floating-point)
  SUBDATASET_5_NAME=NETCDF:"G:\Newman_AK_Data\alaska_daily_1981_050.nc4":pcp
  SUBDATASET_5_DESC=[365x925x1280] pcp (32-bit floating-point)
  SUBDATASET_6_NAME=NETCDF:"G:\Newman_AK_Data\alaska_daily_1981_050.nc4":snow_frac
  SUBDATASET_6_DESC=[365x925x1280] snow_frac (64-bit floating-point)
  SUBDATASET_7_NAME=NETCDF:"G:\Newman_AK_Data\alaska_daily_1981_050.nc4":t_mean
  SUBDATASET_7_DESC=[365x925x1280] t_mean (32-bit floating-point)
  SUBDATASET_8_NAME=NETCDF:"G:\Newman_AK_Data\alaska_daily_1981_050.nc4":t_range
  SUBDATASET_8_DESC=[365x925x1280] t_range (32-bit floating-point)
  SUBDATASET_9_NAME=NETCDF:"G:\Newman_AK_Data\alaska_daily_1981_050.nc4":wetPert
  SUBDATASET_9_DESC=[925x1280] wetPert (64-bit floating-point)
  SUBDATASET_10_NAME=NETCDF:"G:\Newman_AK_Data\alaska_daily_1981_050.nc4":windPert
  SUBDATASET_10_DESC=[925x1280] windPert (64-bit floating-point)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

Note it says the size is 512,512

I made this script to export the variable t_mean in a properly georeferenced format:

    ds = xarray.open_dataset(files[1])
    
    data = ds.t_mean.load().values
    data = data[:,::-1,:]
    data[np.isnan(data)] = -999.0
    time = ds.time.values
    xcoord = pr_targ.x.values
    ycoord = pr_targ.y.values
    
    ds_new = xarray.Dataset({
        'temp': xarray.DataArray(
                    data   = data,   # enter data here
                    dims   = ['time','y','x'],
                    coords = {'time': time,
                              'y':ycoord,
                              'x':xcoord,
                              'spatial_ref':pr_targ.spatial_ref},
                    attrs  = {
                        '_FillValue': -999.0,
                        'units'     : 'deg_C'})},
        attrs= {'long_name' : 'estimated daily mean temperature'})
    ds_new = ds_new.rio.write_crs(pr_targ.rio.crs)
    
    encoding = {"temp": {'zlib': True,'complevel': 9}}
    ds_new.to_netcdf(r"c:\temp\_t_mean_1981_v3.nc", encoding = encoding)

ds_new looks like:

<xarray.Dataset>
Dimensions:      (time: 365, x: 1280, y: 925)
Coordinates:
  * time         (time) datetime64[ns] 1981-01-01 1981-01-02 ... 1981-12-31
  * y            (y) float64 2.37e+06 2.368e+06 ... 5.244e+05 5.224e+05
  * x            (x) float64 -7.342e+05 -7.322e+05 ... 1.822e+06 1.824e+06
    spatial_ref  int32 0
Data variables:
    temp         (time, y, x) float32 -999.0 -999.0 -999.0 ... -999.0 -999.0
Attributes:
    long_name:     estimated daily mean temperature
    grid_mapping:  spatial_ref

When I run gdalinfo on my exported file I get a printout saying there are 365 bands with a size 1280,925. I don't know where all the metadata came from.

C:\Users\Public\Desktop\OSGeo4W>gdalinfo "C:\temp\_t_mean_1981_v3.nc"
Warning 1: Recode from UTF-8 to CP_ACP failed with the error: "Invalid argument".
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Driver: netCDF/Network Common Data Format
Files: C:\temp\_t_mean_1981_v3.nc
Size is 1280, 925
Coordinate System is:
BOUNDCRS[
    SOURCECRS[
        PROJCRS["NAD83 / Alaska Albers",
            BASEGEOGCRS["NAD83",
                DATUM["North American Datum 1983",
                    ELLIPSOID["GRS 1980",6378137,298.257222101,
                        LENGTHUNIT["metre",1]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433]],
                ID["EPSG",4269]],
            CONVERSION["unnamed",
                METHOD["Albers Equal Area",
                    ID["EPSG",9822]],
                PARAMETER["Latitude of 1st standard parallel",55,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8823]],
                PARAMETER["Latitude of 2nd standard parallel",65,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8824]],
                PARAMETER["Latitude of false origin",50,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8821]],
                PARAMETER["Longitude of false origin",-154,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8822]],
                PARAMETER["Easting at false origin",0,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8826]],
                PARAMETER["Northing at false origin",0,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8827]]],
            CS[Cartesian,2],
                AXIS["x",east,
                    ORDER[1],
                    LENGTHUNIT["metre",1]],
                AXIS["y",north,
                    ORDER[2],
                    LENGTHUNIT["metre",1]],
            ID["EPSG",3338]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["Transformation from NAD83 to WGS84",
        METHOD["Position Vector transformation (geog2D domain)",
            ID["EPSG",9606]],
        PARAMETER["X-axis translation",0,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",0,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",0,
            ID["EPSG",8607]],
        PARAMETER["X-axis rotation",0,
            ID["EPSG",8608]],
        PARAMETER["Y-axis rotation",0,
            ID["EPSG",8609]],
        PARAMETER["Z-axis rotation",0,
            ID["EPSG",8610]],
        PARAMETER["Scale difference",1,
            ID["EPSG",8611]]]]
Data axis to CRS axis mapping: 1,2
Origin = (-735223.205799999996088,2371412.932599999941885)
Pixel Size = (2000.000000000000000,-2000.000000000000000)
Metadata:
  NC_GLOBAL#grid_mapping=spatial_ref
  NC_GLOBAL#long_name=estimated daily mean temperature
  NETCDF_DIM_EXTRA={time}
  NETCDF_DIM_time_DEF={365,10}
  NETCDF_DIM_time_VALUES={0,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,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306,307,308,309,310,311,312,313,314,315,316,317,318,319,320,321,322,323,324,325,326,327,328,329,330,331,332,333,334,335,336,337,338,339,340,341,342,343,344,345,346,347,348,349,350,351,352,353,354,355,356,357,358,359,360,361,362,363,364}
  spatial_ref#crs_wkt=PROJCS["NAD83 / Alaska Albers",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",55],PARAMETER["standard_parallel_2",65],PARAMETER["latitude_of_center",50],PARAMETER["longitude_of_center",-154],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","3338"]]
  spatial_ref#false_easting=0
  spatial_ref#false_northing=0
  spatial_ref#geographic_crs_name=NAD83
  spatial_ref#GeoTransform=-735223.2058 2000.0 0.0 2371412.9326 0.0 -2000.0
  spatial_ref#grid_mapping_name=albers_conical_equal_area
  spatial_ref#horizontal_datum_name=North American Datum 1983
  spatial_ref#inverse_flattening=298.257222101
  spatial_ref#latitude_of_projection_origin=50
  spatial_ref#longitude_of_central_meridian=-154
  spatial_ref#longitude_of_prime_meridian=0
  spatial_ref#prime_meridian_name=Greenwich
  spatial_ref#projected_crs_name=NAD83 / Alaska Albers
  spatial_ref#reference_ellipsoid_name=GRS 1980
  spatial_ref#semi_major_axis=6378137
  spatial_ref#semi_minor_axis=6356752.314140356
  spatial_ref#spatial_ref=PROJCS["NAD83 / Alaska Albers",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",55],PARAMETER["standard_parallel_2",65],PARAMETER["latitude_of_center",50],PARAMETER["longitude_of_center",-154],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","3338"]]
  spatial_ref#standard_parallel={55,65}
  spatial_ref#towgs84={0,0,0,0,0,0,0}
  temp#coordinates=spatial_ref
  temp#grid_mapping=spatial_ref
  temp#units=deg_C
  temp#_FillValue=-999
  time#calendar=proleptic_gregorian
  time#units=days since 1981-01-01 00:00:00
  x#_FillValue=nan
  y#_FillValue=nan
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"C:\temp\_t_mean_1981_v3.nc":temp
  SUBDATASET_1_DESC=[365x925x1280] temp (32-bit floating-point)
Corner Coordinates:
Upper Left  ( -735223.206, 2371412.933) (173d29'31.42"W, 70d20'45.14"N)
Lower Left  ( -735223.206,  521412.933) (165d17'28.04"W, 54d 8'34.43"N)
Upper Right ( 1824776.794, 2371412.933) (111d15'32.20"W, 65d49'33.33"N)
Lower Right ( 1824776.794,  521412.933) (127d15' 0.40"W, 51d20'53.50"N)
Center      (  544776.794, 1446412.933) (143d19' 4.19"W, 62d35'42.16"N)

Why are the is the source file gdalinfo size (512,512) and subseted data file (1280,925) different? I must be doing something wrong.

Thanks!

[snowman2]

Note the subdataset mentioned in the original netCDF file:

  SUBDATASET_7_NAME=NETCDF:"G:\Newman_AK_Data\alaska_daily_1981_050.nc4":t_mean
  SUBDATASET_7_DESC=[365x925x1280] t_mean (32-bit floating-point)

It looks like it has the expected dimensions.

You could try opening it up with rioxarray.open_rasterio to get a properly geo-referenced dataset:

import rioxarray

xds = rioxarray.open_rasterio("G:\Newman_AK_Data\alaska_daily_1981_050.nc4", variable="t_mean")