1

I have a CSV file with climate data, latitude, and longitude coordinates.

RH      PREC    TEMP      WS      lat     lon
65.93   1.51    287.43    4.67    63.31   177.62
70.43   1.89    285.32    1.54    63.55   177.98

When I load it into QGIS I get a point shapefile (as expected).

qgis

Right now, everthing is more or less 1D. I know the shape should be something like (76,109), see image above. And then I would add my different data variables to these in the NetCDF file with a resolution of 0.3.

I tried this but got the same scattered plot but now in a NetCDF file. I actually want the points to be rasterised at a 0.3 degree resolution with all the climate data included and then made into a NetCDF file, like this.

NARR resolution

So, how do I convert the points (first image) into a NetCDF file with gridcells of 0.3 degree? I am flexible with program to use, as long as it is not ArcGIS and Java.

1:

8
  • Do you know what grid/projection your data is originally defined on? Neither the input or output samples look like a 0.3x0.3 degree grid. The NCEP NARR Reanalysis grid you show as an example is a 32km Lambert Conformal Conic grid psl.noaa.gov/data/gridded/data.narr.html and "is approximately 0.3 degrees (32km) resolution at the lowest latitude."
    – Dave X
    Commented Nov 29, 2020 at 14:29
  • @Vince, the points are based on NARR, so the centroids should somehow be 0.3 degree separate. The projection is ESPG 102002 Albers Canada conic.
    – Thomas
    Commented Nov 29, 2020 at 18:26
  • what projection is your qgis map using? the coordinates at the bottom make me think this isn't lat-long. if you change the project CRS (not the layer CRS) to epsg:4326 (lat-long) then do you see a regular, lat-long aligned grid? The layer CRS, judging by your CSV snippet, should be set by you to epsg:4326 when you import the CSV.
    – Spacedman
    Commented Nov 30, 2020 at 10:52
  • if the points aren't on the grid system you want then you have two problems: 1. imterpolating your data to the grid system you want and 2. saving that regridded data to a NetCDF.
    – Spacedman
    Commented Nov 30, 2020 at 10:55
  • 1
    What do you mean by "buffered with the exact 0.3 degree"? You've got a set of irregular points in lat-long that might not be anywhere near a regular 0.3 degree grid. To transpose the values to a regular 0.3 degree grid involves techniques like averaging the values at points near to each 0.3 degree grid point, aka "spatial interpolation". Is that what you need to do?
    – Spacedman
    Commented Nov 30, 2020 at 11:26

1 Answer 1

1

The NARR grid isn't a 0.3 degree grid, is defined directly as a regular 32.463km regular grid with square pixels in a projected Lambert Conic Conformal projection space. You should work with your data in its native projection.

You can open up a NARR GRIB file downloaded from https://www.ncei.noaa.gov/thredds/fileServer/model-narr-a-files/201410/20141002/narr-a_221_20141002_0000_000.grb and with gdalinfo get it's raster georeferencing info:

$ gdalinfo ~/Downloads/narr-a_221_20141002_0000_000.grb |head -60
Driver: GRIB/GRIdded Binary (.grb, .grb2)
Files: /Users/drf/Downloads/narr-a_221_20141002_0000_000.grb
Size is 349, 277
Coordinate System is:
PROJCRS["unnamed",
    BASEGEOGCRS["Coordinate System imported from GRIB file",
        DATUM["unnamed",
            ELLIPSOID["Sphere",6371200,0,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["unnamed",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",50,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",-107,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",50,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",50,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        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["easting",east,
            ORDER[1],
            LENGTHUNIT["Metre",1]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["Metre",1]]]
Data axis to CRS axis mapping: 1,2
Origin = (-5648873.725474949926138,4363473.848625719547272)
Pixel Size = (32463.000000000000000,-32463.000000000000000)
Corner Coordinates:
Upper Left  (-5648873.725, 4363473.849) (148d23'40.62"E, 46d30'58.42"N)
Lower Left  (-5648873.725,-4628777.151) (145d32'24.16"W,  0d51'19.97"N)
Upper Right ( 5680713.275, 4363473.849) (  2d19'25.78"W, 46d14' 4.78"N)
Lower Right ( 5680713.275,-4628777.151) ( 68d16'45.97"W,  0d45'10.94"N)
Center      (   15919.775, -132651.651) (106d46'57.60"W, 48d48'25.02"N)

...

Or you can open it directly in QGIS and get it's projection info as +proj=lcc +lat_0=50 +lon_0=-107 +lat_1=50 +lat_2=50 +x_0=0 +y_0=0 +R=6371200 +units=m +no_defs +type=crs

If you want to create a NARR-like NetCDF file, you could start with a template file from https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/north-american-regional-reanalysis-narr like (drill down trhough the TDS server to a file, and then choose the NetCDFsubset option to get: https://www.ncei.noaa.gov/thredds/catalog/model-narr-a-files/201410/20141002/catalog.html?dataset=model-narr-a-files/201410/20141002/narr-a_221_20141002_0000_000.grb and choose a couple compatible variables. (I'd choose a couple surface variables, disable the horizontal subsetting, and enable the CF-compliant lat/lon vars. I'd also choose a couple simple surface-level variables as templates.) Then use python with NetCDF4 to update and replace the data in those variables and their metadata with your own.

NetCDF subset selection screenshot

Alternately, you could use the projection information above to reformat your data into a stack of ArcInfo ASCII Grids .AAG files, and then use GDAL to translate them into a multiband TIFF or stick-built NetCDF files with something like:

gdal_merge -o myfile.tif -of GTiff -separate RH.aag PREC.aag TEMP.aag  WS.aag 

gdal_merge -o myfile.nc -of netCDF -separate RH.aag PREC.aag TEMP.aag  WS.aag 

Oh, you'd also need RH.prj, PRED.prj, TEMP.prj, and WS.prj files with the WKT of the data:

PROJCS["unnamed",GEOGCS["GCS_Coordinate_System_imported_from_GRIB_file",DATUM["D_unnamed",SPHEROID["Sphere",6371200.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-107.0],PARAMETER["Standard_Parallel_1",50.0],PARAMETER["Standard_Parallel_2",50.0],PARAMETER["Latitude_Of_Origin",50.0],UNIT["Meter",1.0]]

For the full-size grb file data, the .aag files https://en.wikipedia.org/wiki/Esri_grid#ASCII would have headers like this from:

gdal_translate  -of AAIGrid -b 1 ~/Downloads/narr-a_221_20141002_0000_000.grb MSLET_Pa.aa

ncols        349
nrows        277
xllcorner    -5648873.725474949926
yllcorner    -4628777.151374280453
cellsize     32463.000000000000
NODATA_value  9999
 9999.0 9999 9999 9999 9999 9999 9999 101239.625 101178.625 101113.625 101049.625 100987.625 100925.625 100861.625 100796.625 100728.625 100659.625 100594.625 100534.625 100479.625 100432.625 100392.625 100360.625 100333.625 100308.625 100285.625 100263.625 100238.625 100209.625 100178.625 100148.625 100115.625 100076.625 100034.625 99991.625 99951.625 99912.625 99877.625 99846.625 99833.625 99837.625 99857.625 99898.625 99956.625 100015.625 100094.625 100192.625 100292.625 100391.625 100497.625 100604.625 100710.625 100809.625 100898.625 101005.625 101110.625 101198.625 101285.625 101375.625 101459.625 101536.625 101601.625 101660.625 101724.625 101789.625 101844.625 101906.625 101988.625 102082.625 102184.625 102259.625 102322.625 102382.625 102436.625 102482.625 102526.625 102562.625 102596.625 102626.625 102648.625 102665.625 102679.625 102690.625 102697.625 102698.625 102689.625 102671.625 102646.625 102611.625 102572.625 102526.625 102474.625 102420.625 102365.625 102310.625 102255.625 102201.625 102146.625 102095.625 102044.625 101993.625 101945.625 101896.625 101848.625 101802.625 101757.625 101711.625 101666.625 101620.625 101576.625 101534.625 101492.625 101453.625 101417.625 101383.625 101352.625 101325.625 101301.625 101280.625 101262.625 101248.625 101236.625 101225.625 101216.625 101208.625 101201.625 101193.625 101188.625 101182.625 101177.625 101172.625 101167.625 101162.625 101158.625 101153.625 101150.625 101148.625 101146.625 101145.625 101146.625 101148.625 101152.625 101156.625 101161.625 101167.625 101175.625 101183.625 101193.625 101204.625 101216.625 101230.625 101245.625 101261.625 101279.625 101297.625 101315.625 101335.625 101356.625 101377.625 101398.625 101420.625 101441.625 101461.625 101482.625 101501.625 101520.625 101537.625 101555.625 101572.625 101588.625 101603.625 101617.625 101631.625 101644.625 101656.625 101667.625 101677.625 101687.625 101696.625 101705.625 101713.625 101720.625 101727.625 101733.625 101739.625 101744.625 101748.625 101751.625 101754.625 101757.625 101759.625 101761.625 101762.625 101763.625 101762.625 101762.625 101760.625 101758.625 101755.625 101753.625 101750.625 101747.625 101741.625 101733.625 101726.625 101714.625 101716.625 101723.625 101728.625 101725.625 101728.625 101731.625 101731.625 101734.625 101728.625 101713.625 101685.625 101643.625 101602.625 101559.625 101511.625 101462.625 101417.625 101369.625 101314.625 101245.625 101168.625 101088.625 101007.625 100929.625 100867.625 100800.625 100735.625 100677.625 100614.625 100529.625 100475.625 100439.625 100398.625 100337.625 100239.625 100152.625 100081.625 100012.625 99951.625 99916.625 99871.625 99815.625 99741.625 99655.625 99568.625 99479.625 99392.625 99311.625 99243.625 99179.625 99114.625 99049.625 98984.625 98922.625 98865.625 98814.625 98769.625 98733.625 98706.625 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
...

You could fill your data into a portion of the whole array pre-filled with nodata, or you could figure out the subset of the data you do have and change the header to match.

7
  • Thanks @Dave X. When I try to get the data from unidata as you provided. I only get the following message Validation errors: {thredds.server.ncSubset.validation.start_gt_end} -- . And, I cannot seem to subset for the lat lon I need. Do you know how to fix it?
    – Thomas
    Commented Nov 30, 2020 at 19:45
  • What area do you need? The lat/lon ranges can be sort of odd in a projected system like this. A quick look at the NARR domain in QGIS has the data crossing the dateline, and the middle of the grid 40 degrees higher than the corners.
    – Dave X
    Commented Nov 30, 2020 at 20:02
  • I need data for Alaska, Yukon & Northwest Territories, give me a sec for the lat lon, please.
    – Thomas
    Commented Nov 30, 2020 at 20:07
  • I grabbed the full file. The start_gt_end bit might be interacting with the 180W dateline. Maybe a bounding box of 220W might work.
    – Dave X
    Commented Nov 30, 2020 at 20:13
  • The bounding box would be something like: (180, 85) (259, 85) (180, 48) (259, 48)
    – Thomas
    Commented Nov 30, 2020 at 20:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.