Export netCDF data based on latitude and longitude

Hello,

I have the following code which coded by KSSV This code is to read netCDF and export them to excel file. The problems are:

1. I want to export a specific long, lat
2. export the data to several csv files instead one excel file

I appreciate your help and time!

files = dir('*.nc') ;

nfiles = length(files) ;
P = cell(nfiles,1) ; % precipitation of all files
date = cell(nfiles,2) ;
time = cell(nfiles,2) ;
for i = 1:nfiles
    %%Read dat
    date{i,1} = ncreadatt(files(i).name,'/','HDF5_GLOBAL.BeginDate')  ;
    date{i,2} = ncreadatt(files(i).name,'/','HDF5_GLOBAL.EndDate')  ;
    %%Read time
    time{i,1} = ncreadatt(files(i).name,'/','HDF5_GLOBAL.BeginTime')  ;
    time{i,2} = ncreadatt(files(i).name,'/','HDF5_GLOBAL.EndTime')  ;
    %%REad precipitation
    P{i} = ncread(files(i).name,'precipitation') ;
end

myfile = 'myfile.xlsx' ;
for i = 1:nfiles
    xlswrite(myfile,P{i},i)
end

Answer 1

[Walter Roberson]

You can use the same kind of masking techniques I talk about in https://www.mathworks.com/matlabcentral/answers/347505-extracting-netcdf-data-based-on-time#answer_272996

You have the advantage of not having to make a correction for the odd way netcdf encodes time.

[Majid Mohamod]

I appreciate your answer. I tried this script and it give me an error "it's included with this question"

clear all
close all
ncid = netcdf.open('uwnd.sig995.2016.nc','NC_NOWRITE');
lat = netcdf.getVar(ncid,0);
whos lat
lon = netcdf.getVar(ncid,1);
whos lon
time = netcdf.getVar(ncid,2);
whos time
wind = netcdf.getVar(ncid,3); %var3 wind 4x day
whos wind
latlim = [5 -5]; %5N 5S
lonlim = [-50 -45]; %-50W -45W

I like to share some points with you: 1- my files are already separated based on time "each nc4 file is one day" 2- My nc4 files cover the whole world. 3- I want to read all my nc4 file "which is done so far", extract my study area data and output each nc4 to Excel file as Lat, Long and Precipitation value.

What I have achieved so far, that I can read the whole nc4 files and export the global data as spreadsheets in one excel file which become difficult to open it due to the size of excel file.

This is the script what I used:

   files = dir('*.nc4') ;

nfiles = length(files) ;
P = cell(nfiles,1) ; % precipitation of all files
date = cell(nfiles,2) ;
time = cell(nfiles,2) ;
for i = 1:nfiles
    %%Read dat
    date{i,1} = ncreadatt(files(i).name,'/','HDF5_GLOBAL.BeginDate')  ;
    date{i,2} = ncreadatt(files(i).name,'/','HDF5_GLOBAL.EndDate')  ;
    %%Read time
    time{i,1} = ncreadatt(files(i).name,'/','HDF5_GLOBAL.BeginTime')  ;
    time{i,2} = ncreadatt(files(i).name,'/','HDF5_GLOBAL.EndTime')  ;
    % REad precipitation
    P{i} = ncread(files(i).name,'precipitation') ;
end

myfile = 'myfile.xlsx' ;
for i = 1:nfiles
    xlswrite(myfile,P{i},i)
end

Do you think should I merge the nc4 files first by using this script, if so, I tried but it give me an error too!

the script is:

files = dir('*.nc4') ;

nfiles = length(files) ;
P = cell(nfiles,1) ; % precipitation of all files
date = cell(nfiles,2) ;
time = cell(nfiles,2) ;
for i = 1:nfiles
    %%Read dat
    date{i,1} = ncreadatt(files(i).name,'/','HDF5_GLOBAL.BeginDate')  ;
    date{i,2} = ncreadatt(files(i).name,'/','HDF5_GLOBAL.EndDate')  ;
    %%Read time
    time{i,1} = ncreadatt(files(i).name,'/','HDF5_GLOBAL.BeginTime')  ;
    time{i,2} = ncreadatt(files(i).name,'/','HDF5_GLOBAL.EndTime')  ;
    %%REad precipitation
    P{i} = ncread(files(i).name,'precipitation') ;
end

%%Write all the above data into nc-file
ncfile = 'merged.nc' ;
[nx,ny] = size(P{1}) ;
%
nccreate(ncfile,'time','Dimensions',{'time',1,Inf},'DeflateLevel',8) ;
%
nccreate(ncfile,'precipitation','Dimensions',{'lat',nx,'lon',ny,'time'},'DeflateLevel',8) ;
ncwriteatt(ncfile,'precipitation','units','mm');
ncwriteatt(ncfile,'precipitation','long_name','Daily accumulated precipitation (combined microwave-IR) estimate with gauge calibration over land');
ncwriteatt(ncfile,'precipitation','coordinates','lat lon');
ncwriteatt(ncfile,'precipitation','FillValue','-9999.9004');
ncwriteatt(ncfile,'precipitation','origname','precipitation');
ncwriteatt(ncfile,'precipitation','fullnamepath','/precipitation');
for i = 1:nfiles
    ncwrite(ncfile,'precipitation',P{i},[1,1,i]) ;
end
ncwriteatt(ncfile,'/','HDF5_GLOBAL.BeginDate',date{1,1}) ;
ncwriteatt(ncfile,'/','HDF5_GLOBAL.BeginTime',time{1,1}) ;
ncwriteatt(ncfile,'/','HDF5_GLOBAL.EndDate',date{nfiles,2}) ;
ncwriteatt(ncfile,'/','HDF5_GLOBAL.EndTime',time{1,2}) ;
ncwriteatt(ncfile,'/','creation_date',datestr(now));
ncdisp(ncfile

Your script is so helpful if i have combined nc4 files, but I don't have such file.

   from_year = 2014; to_year = 2015;   %for example

        time = netcdf.getVar(ncid,2);
        wind = netcdf.getVar(ncid,3); %var3 wind 4x day

        time_datenum = time / 24 + datenum('1800-01-01 00:00:0');
        dv = datevec(time_datenum);
        date_match = dv(:,1) >= from_year & dv(:,1) <= to_year;

        selected_wind = wind(:, :, date_match);

        Or, expanding the selection

        from_date = '2014-04-01'; to_date = '2015-03-31';

        time = netcdf.getVar(ncid,2);
        wind = netcdf.getVar(ncid,3); %var3 wind 4x day

        time_datenum = time / 24 + datenum('1800-01-01 00:00:0');
        date_match = time_datenum >= datenum(from_date) && time_datenum <= datenum(to_date);

        selected_wind = wind(:, :, date_match);

These all I've been stuck in this problem. I will be so grateful to solve this problem. I really need these data in my research.

[Walter Roberson]

I see two difficulties with your file:

1) There is no clear mapping between latitude and longitude and the dimensions (480 x 1440). Are we to assume that the latitudes are 0.375 degree apart and the longitudes are 0.25 degree apart?

2) There is no geo-reference to handle projection from a rectangular array to a (more or less) spherical array.

We could make some guesses to handle problem #1, but it is not immediately obvious that you could safely take a rectangular subset of the data and still produce a "fair" precipitation table.

[Majid Mohamod]

Please, can we just assume the latitudes are 0.375 degree and longitudes are 0.25 degree. I can fix it later once I get access to download the new data, "I am working on that trying to figure out why my data generate in the way".

For the projection, can you explain more, you mean whenever we convert projection, we will get an error in the accuracy. If so, I don't think it will be matter if the error is small.

I also attached some information about the these data product. I hope it would help.

 Thank you a lot

[Walter Roberson]

Do you have the Mapping Toolbox?

[Majid Mohamod]

I don't have it "Mapping Toolbox" I am going to search about it.

I want to keep you updated, there are some helpful information below:

This is what I found about the Horizontal Resolution:

[Walter Roberson]

At 50N, the width of a degree of longitude is less than 45 miles, compared to more than 69 miles at the equator. That is roughly 2/3, which is enough to make a fair difference in analyzing precipitation in any absolute sense.

For example I happen to live at 49.97 N, and the area I live in is just slightly above the cutoff for being a semi desert (there is some desert about 80km west of me). We get rain in the summer but our winter is so cold that it is 5 times dryer than the Majove. If our readings were underestimated by 1/3 then we would appear to be true dessert, which is not at all the case. If our readings were overestimated by 1/3 then we would appear to be substantially wetter than we are. So the correction does matter.

[Majid Mohamod]

I completely agree with you. Any suggestion how we can build this script. I really appreciate your experience.

[Walter Roberson]

You would like to extract a subsection based on latitude and longitude, but what would you like to do with it after that? Do you need to regrid into squares of equal area?

[Walter Roberson]

I have asked someone for some regridding advice.

[Majid Mohamod]

I want to combine them in sequence data. So I will have satellite weather stations. Each of them is composed from (Lat, Long, Precipitation value) for period of time (e.g. a month or a year). If we combine them all (the Global data) to have daily or monthly precipitation (e.g. rainfall). The excel file size will be huge so it will be difficult to work on them. So, it will be better to cut only the interested area from the global scale data, and combine them!

[Walter Roberson]

Unfortunately my contact indicates that they will have a bad internet connection for a few days.

The nc4 files that you posted do not have the position headers that you document above. That would be okay if you could be sure that all of them were the same position.

With the information you have provided, Yes, we can select the data by lat and long range, and write it to an Excel file. I am concerned about how you process the data after that. The grid locations will not be of equal area. As long as you only compare each grid location to itself (e.g., is it getting wetter or drier compared to itself) then that is acceptable. If, however, you need to compare grid locations to each other then you might need to correct for area. The data in its existing form is not suitable for directly computing "rainfall per square kilometer" for example.

[Majid Mohamod]

If we can just read part of the precipitation array even based on the data position in the array (nrow, ncol) I can convert them later to coordinate. I hope if possible to output the data into different excel files instead into one excel file within multiple sheets. I will do some validations on the data later to make sure the coordinate is right.

I appreciate your help!

[Walter Roberson]

%setup done once
lat_vec = linspace(-50, 50, 400);
long_vec = linspace(-180, 180, 1440);
lat_res = 0.25;
long_res = 0.25;

latlim = [5 -5]; %5N 5S
lonlim = [-50 -45]; %-50W -45W

matching_lat = lat_vec > min(latlim) - lat_res/2 & lat_vec < max(latlim) + lat_res/2;
matching_long = long_vec > min(lonlim) - long_res/2 & long_vec < max(lonlim) + long_res/2;

%and for each file,

data_subset = P{i}(matching_lat, matching_long);

The lat_res and long_res are there to take into account that the coordinates are the centers of the grid cells, so the given limits might happen to be within the logical confines of a grid cell and yet not happen to match based upon the given limit. For example, if the limit given were -5.1 then you need to include the grid entry that is nominally at -5 because that grid entry really runs from -5.25 to -4.75 (I do not know whether that is semi-open or semi-closed)

[Majid Mohamod]

Thank you so much man. I will try it and keep you updated. :)