Here's some python code that does what you want, reading GDAL files that represent data at specific times and writing to a single NetCDF file that is CF-Compliant
Convert a bunch of GDAL readable grids to a NetCDF Time Series.
Here we read a bunch of files that have names like:
These kind of files are bit of a pain, because they are using one set of variables to index another set of variables. But once you've figured out what is indexing what, you can do it pretty easily using python. Here we use netCDF4 to read the netcdf and Pandas to write the CSV, but you could use other tools to do the job.
from pylab import *
To get the appropriate cell coordinates from your latitude and longitude you need to know the coverage of the netCDF file (for example, the coordinate of the upper left cell and the size of the cells in the x and y directions). @Luke 's solution below works for a srtaight x/y raster with no cell rotation. If you don't have that, you could look at the affine ...
I figured out the answer some time ago but I hope this will help other people that are dipping their toes into the NetCDF file format.
I found that the easiest way to create a new NetCDF file containing several GeoTIFF files is using GDAL and then using NCO and CDO tools (Provided by unidata) to add the appropriate metadata.
First you convert the GeoTIFF´s ...
Xarray is (intentionally) ignorant of coordinate systems, so it has no special handling for cyclic coordinates such as longitude. Because your longitude array has only increasing values, xarray interprets selections like slice(40, -80) in the same way that x[i:j] works if x is a NumPy array and i > j >= 0, and thus returns an empty selection.
If you want to subset some data from a NetCDF file using a lon/lat bounding box and that NetCDF file is not aligned with east/north, one strategy is to use a point-in-polygon routine and then find the min/max i,j indices of those points to define a subset to extract.
There are several different python packages to deal specifically with netCDF data from ...
Here's what we've been doing to set up THREDDS Data Server (TDS) catalogs for regional oceanographic modeling providers in the US Integrated Ocean Observing System to serve their models results.
There are four basic types of catalogs we have been setting up:
A top level catalog that points to other catalogs that you want exposed
An "all" catalog that ...
The Norwegian Met office has a THREDDS server at http://thredds.met.no/thredds/
so if you see the forecast you are trying to access there, you can extract just the subset you want from the OPeNDAP URL, which NetCDF4-Python treats like a local netcdf file.
url = 'http://thredds.met.no/thredds/dodsC/arome25/...
The metadata of the netcdf file includes projection information, which can be identified by QGIS as a custom CRS:
+proj=lcc +lat_1=50 +lat_2=50 +lat_0=50 +lon_0=-107 +x_0=5632642.22547 +y_0=4612545.65137 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
but there is no extent given in projection units, only latlon for the corners.
If you run gdalinfo ...
So you have 3D data (lon,lat,time) in a netcdf file and you want to extract a time-series as a specific location in Python, right?
Here's one way using netCDF4:
But it's way easier using Xarray:
import xarray as xr
Use this to get access to all bands
lradin <- brick("lrad_ITPCAS-CMFD_V0106_B-01_197912.nc")
Use ?subset or lradin[] to get specific band/s.
Use ?getZ to find out the specific date for each band. Use extract(lradin, cbind(lon, lat)) to get time series at specific point/s.
Feel free to ask other focussed questions like this as needed.
There are ...
Yes, the CF 1.6 Conventions for NetCDF include the specification of collections of time series and it seems your data is similar to example H.2.1 "Orthogonal multidimensional array representation of time series":
If you store your data this way, IDV should be able to recognize this as "point data". Hopefully more applications in the future will take ...
You can shift the data into the correct position using these GDAL commands:
gdal_translate -a_srs EPSG:4326 -a_ullr 0 0 360 -90 NETCDF:"sic_average_nclimate.nc":SIC_Change change360.tif
gdalwarp -t_srs WGS84 change360.tif change180.tif -wo SOURCE_EXTRA=1000 --config CENTER_LONG 0
(with a little help by Frank Warmerdam: How to reproject raster from 0 360 ...
I'd suggest trying NCO again, as it's well suited for this type of operation, and can be done with one command. The tool to concatenate several files along the record dimension (e.g. time) is ncrcat. For example, to combine all the files for May 2018:
ncrcat gdaps_anal_201805????.nc gdaps_anal_201805.nc
The resulting file gdaps_anal_201805.nc should ...
I would recommend looking into rioxarray for your dataset.
You can open your dataset like so:
xds = xarray.open_dataset('D:\Weather data\et_01012018.nc')
If your CRS is not discovered, you should be able to add it like so:
Then, you should be able to create a geotiff ...
Here is an example of using rioxarray to mask out data with a shapefile:
from shapely.geometry import mapping
MSWEP_monthly2 = xarray.open_dataarray('D:\G3P\DATA\Models\MSWEP\MSWEP_monthly.nc4')
These data are not in a netcdf format but, rather a space delimited ASCII format. This data is a bit difficult to deal with because of the lack of headers and any type of unique station identifier.
Returns a vector of file names on disk and pulls associate years from file name.
f <- list.files(getwd())
y <- as.numeric(unlist(lapply(strsplit(f, "[.]")...
I think the issue here might be the way the data is imported not the projection.
If you use the Multidimension Tools tool box in ArcGIS 10.0, you need to enter the dimensions of your netcdf variable as the "Row Dimension" option in order to plot the value at each coordinate.
A big question is whether you are going to read the entire raster from the file into memory before processing it, or whether the file is so large that you will process it incrementally, or process some subset of the overall file.
If you will load it all into memory, then you will be doing mostly sequential access, and the fastest format will be a tossup ...
I know this is an old question, but if anyone is still looking for an answer on how to do this in R the following code that will make a netCDF with 2 variables (prec, tmax), and 3 dimensions (long, lat, time). Pretty straight forward to change this suit your number of variables or time dimension etc.
prec <- getData("...
Thanks to gene mentioning the netcdfbrowser plugin (see the addLayer()-method in netcdfbrowserdialog.py) I managed to get it working. In case of NetCDF-files the file name is not sufficient to load a raster. We also need to include the variable of interest. The following works:
Yes, the NetCDF CF Metadata Conventions version 1.6 specifies how to store point and station time series data in chapter 9 "Discrete Sampling Geometries".
Since your data has the same sample times for all stations, I agree with Rich that you can base your netCDF structure on the example in section H.2.1 "Orthogonal multidimensional array representation of ...
If you have moved on to the netCDF Java API it might be interesting for you to try Python as well. I have worked with Java for a while but find Python for number crunching and data analysis more convenient.
Geoserver should be able to serve your postGIS data as NetCDF. You can then use Python or Java to generate a Netcdf file that contains the difference in ...
I have found the answer to my own question. If this works for anybody else who is having the same problem, maybe they can vote up the question or the answer, since I cannot get any points for answering my own question. :)
When I was looking for an answer I saw that populating the "row dimension " in ArcGIS with the dimensions of the netcdf ...
It might be easier to see what is going on when you create a RasterBrick (multiple layers), either representing depth, or time. The number of layers should make it clear which is is. the z-values too.
You could start with
b <- brick("data.nc", varname="nitrat")
and go from there.
The error in hour ncdf code is because you mix up the ncdf package (open....
I guess that my answer from your other StackOverflow question did not lead you in the right direction?
Here is a more detailed answer that may nudge you in the right direction.
First, we need to know the projection of your data and the extent to be able to project to the Long/Lat grid correctly. Unfortunately, we do not have the PROJ4 CRS or the extent, so ...
This is not really a raster in longitude latitude, it's just arrays of values (including longitude and latitude). You can deal with these explicitly like this:
f <- "F18-SSUSI_EDR-NIGHT-DISK_DD.20150107_SN.26920-00_DF.NC"
## treat these not as rasters, but as arrays of values
## though raster() is extremely helpful in simplifying the ...
We have one of these files on hand:
## dp is the root to our local data repository
f <- file.path(dp, "data", "ftp.cdc.noaa.gov/Datasets/ncep.reanalysis2.derived/pressure/air.mon.mean.nc")
b <- brick(f, level = 1)
On my system that uses the "ncdf4" package, but it could also use "ncdf".
The global Pacific-view extent looks close, ...
Your problem comes from adressing the subdatasets wrong.
If you run gdalinfo on the complete file it will display the names of the subdatasets:
To get the information of the first subdataset you need to feed the complete name into gdalinfo