subscript indices must either be real positive integers or logicals - how to change lat/lons from double to integers?

Hi,

I am executing a script for a global grid, using for i=1:144 for j=1:91 [z]=SPI ...etc end end

But I get the error message: subscript indices must either be real positive integers or logicals

The lat and lons were originally 'double' format in the source netcdf file... is there an easy way to change them to integers in matlab?

Answer 1

[Geoff]

If they just need to be indices, you can use ceil, floor, fix or round. eg

x = fix(x)

If you want to convert them to ACTUAL integer types, you can use:

x = uint32(x);

[EDIT] As the question changes...

As Walter pointed out, you are overwriting Z every loop. Perhaps you wanted Z indexed by i and j... In which case:

Z(i,j) = SPI(precip30y_model1(i,j),12,12);

But since you said that Z is empty after iteration, is it possible that SPI() can return an empty set instead of a scalar? In that case:

Z = nan(144,91);
for i=1:size(Z,1)
  for j=1:size(Z,2)
    s = SPI(precip30y_model1(i,j),12,12);
    if ~isempty(s)
      Z(i,j) = s;
    end
  end
end

That uses NaN to represent "no data".

[EDIT 2012-05-23]

Okay, from your comments I think you want this:

Z = nan(144,91);
for i=1:size(Z,1)
  for j=1:size(Z,2)
    Z(i,j) = SPI(squeeze(precip30y_model1(i,j,:)),12,12);
  end
end

[Madeleine P.]

Hi, well that solved the problem (thanks for that!), only it revealed another problem; the script (not mine) I'm implementing is:

[Z]=SPI(precip30yTEMLS(i,j),12,12);

But I'm now getting an error message: Attempted to access SPI(0.166925,12,12); index must be a positive integer or logical...

But it's the script (SPI.m) which is calculating the index (SPI), and it's complaining that its in the wrong format? Any ideas on what I can do to fix this (other than email the author of the script)? (b.t.w I'm new to matlab -normally use ncl...)

[Geoff]

Do you by any chance have a variable in your workspace called SPI? Try:

clear SPI

The error message is caused by your value precip30yTEMLS(i,j) being a non-integer, and apparently you are trying to use this to index a matrix called SPI.

[Madeleine P.]

Hi again, yep- this worked (I had previously created a matrix SPI). So now I can run my script for the globe:

for i=1:144

for j=1:91

[Z]=SPI(precip30y_model1(i,j),12,12);

end

end

.... and the command executes without any error messages, but creates an empty array (?!) i.e. in the workspace window, I've got this array Z[] (i.e. empty, when it should have 144x91 values). The input dataset on which the SPI calculation is performed is 144x91x360 (lon,lat,precipitation). Not sure what I am doing wrong....hard to decipher if I'm not getting any error msgs. Any ideas ? Thanks.

[Walter Roberson]

That loop is the same as calculating *only*

SPI(precip30y_model1(144,91),12,12)

as you are overwriting Z on every trip through the loop.

[Madeleine P.]

I think I know what's wrong, but don't know how to fix: my input dataset is 3 dimensions of rainfall, time (360), lat (91), lon(144) --> when ncread into matlab it gets transformed into lon,lat,time. This is monthly data which gets is supposed to be transformed into 1 output variable (SPI) per gridcell (i,j,) (via the gamfit and gamcdf matlab functions). So I am currently using 3 'for loops' i.e. i=144, j=91 and k=360, (without the k for loop, I get an error message 'subscripted assignment dimension mismatch', but Z cannot be calculated every month - it needs years of data (at least 30) to calculate 1 value of SPI, so if I have the time loop getting the script to try to calculate SPI every month, it's not going to work... this is why the script says the data should be month data vector not matrix (monthly or seasonal precipitation).

How to fix this, given that my input data is in 3 dimensions -I cannot use the time for loop, but get the above error message when I do... ?

[Geoff]

Are you saying that 'precip30y_model1' is the 3D matrix and you want to pass all readings from a single lat/lon into the SPI function as a vector?

In that case you could use the colon operator on the third dimension and squeeze the result... like this:

SPI( squeeze(precip30y_model1(i,j,:)), 12, 12 )

[Madeleine P.]

Well, that does work, if I sandwich it between the lon,lat,time for loops, I am getting Z being calculated, but there is something weird going on -after 7.5 hours of runtime, I had a Z array of 1x348, which doesn't really make sense, as I am trying to create a lon/lat array of 144x91....

I am using [Z]= SPI(squeeze(precipmodel1(i,j,:)),12,12); ...as in the workspace window, it's telling me that this precipmodel variable has 144x91x360 dimensions (so the '360' dimension is the one I want to 'squeeze'). However, it appears that the time(360) variable is not the variable being squeezed and it's either lat or lon... I'll try [Z]= SPI(squeeze(precipmodel1(:,j,k)),12,12); and see if that works....

[Geoff]

The code I gave in my last comment will give you a vector of 360 values for the long / lat at i / j.

You need to assign it to Z(i,j), not Z.

[Geoff]

Have edited the answer. I _think_ that's what you're trying to achieve.

[Madeleine P.]

Hi- I tried your code (above), and immediately got an error message 'subscripted assignment dimension mismatch' ... :-(

[Geoff]

Do you know how to use the debugger? Set a breakpoint on that line (F12), run your code, execution breaks. F10 steps through code line by ine. F5 runs. You can find out a lot that way...

That said, what is the result of this: size(SPI(squeeze(precip30y_model1(1,1,:)),12,12))

I suspect it will be [1,348]. That's what you got when you were overwriting Z every iteration. The question then will be: are you conceptualising your problem incorrectly, is the SPI function broken, or do you need to do something to that (ie mean or sum) result before you put it into your Z array?

Interestingly, 348 is 360 - 12, and you seem to be passing that number into the function a couple of times.

[Madeleine P.]

Ok, some progress has been made, at least in understanding of what I want to calculate. I do actually calculate a value of Z for each lon,lat,time values. So Z - which is a function, should actually be 3 dimensional. So I ran your code above (2012-05-23) for 3 dimensions (and no squeeze) and I got the error message 'subscript indices must either be real positive integers or logicals.

My code was:

Z = nan(144,91);

for i=1:size(Z,1)

for j=1:size(Z,2)

for k=1:size(Z,3)

Z(i,j,k) = SPI(precip30y_model1(:,:,:)),12,12);

end

end

end

When I tried the same thing with just [Z] = SPI(precip30y_model1(:,:,:)),12,12); ... I get the error message 'colon operands must be real scalars, yet when I use the [Z] = SPI(precip30y_model1(i,j,k)),12,12); I get an empty array [Z] for my troubles, as we've already covered previously.

And yet, when I run the script for 1 gridpoint:

[Z] = SPI(precip30y_model1(116,87,:)),12,12);

I get an array of [Z] = 1x348 (which is now what I realise I should get) -and the values all look right (and 348 is ok too as the script discards the first lot of months according to the scale (i.e. in this case 12). So the script works for 1 point but I just can't get it to work for the entire globe.

If I create the NaN array, then use the 3 loops for lon, lat, time, then [Z]=SPI(precip(144,91,360),12,12); I get a blank array really quickly (i.e. Z is not array, and yet I have never made [Z(i,j,k)] work.

So at this point, I'm at a loss. I hear NCL are bringing out their own function for SPI sometime this month - maybe I'll have better luck with it.

Anyway, thanks for your help- I really appreciate it. However I think that this function was not designed to be looped over a lat/lon grid....

[Walter Roberson]

Try this experiment:

S = size(precip30y_model1)

Z = nan(348,S(1),S(2));

for j=1:S(1)

for k=1:S(2)

Z(:,j,k) = SPI(precip30y_model1(j,k,:)),12,12);

end

end

[Geoff]

Right, so it's not a question of having requested the right data. It's that you are trying to store vectors in a single element of a matrix. I agree with Walter's suggestion (expand your Z matrix to 3 dimensions). Alternatively, you could use cell arrays: declare Z = cell(144,91) and use curly braces when assigning: Z{i,j} = SPI( etc ), but that might make life difficult later. In the case of Walter's suggestion, I'd consider putting your results into dimension 3 instead of dimension 1, so that it is consistent with your precip30y_model matrix.

[Walter Roberson]

Putting the results in dimension 1 is much more efficient. You can always permute(Z,[2 3 1]) afterwards

[Madeleine P.]

Hi, this seemed to work, in that I have an array of Z(348,144,91), and when I try to look at the data max(Z) and min(Z) I get reasonable answers. (and not nans). Is there like an ncdump command for matlab which will let me look at all the data stored in Z? Unless I can do this i won't know if I am calculating exactly what I want. I have tried to write it to a netcdf file so I can export it and work with it in a language I'm actually familiar with (i.e. NCL), but am having trouble; I'm trying:

> nccreate('spi_model1.nc','spi');

>ncwrite('spi_model1.nc','spi',Z);

But am getting an error message 'number of input elements does not match variable size'.

I tried the harder way:

spi=Z(360,144,91);

ncid=netcdf.create('spi_model1.nc','NOCLOBBER');

dimid=netcdf.defDim(ncid,'time',348);

dimid=netcdf.defDim(ncid,'lon',144);

dimid=netcdf.defDim(ncid,'lat',91);

varid=netcdf.defVar(ncid,'spi','NC_DOUBLE',dimid);

netcdf.putVar(ncid,varid,spi);

netcdf.close(ncid);

---this got me the same error message! And yet, in my workspace window, it's telling me I have a variable Z(348,144,91) ? It would seem like the variable Z may not contain what it's supposed to contain; I need to look at all the data -how do I do this?

[Walter Roberson]

nccreate() creates a scalar variable unless you specify dimension information.

To fix the other problem:

dimid(1) = netcdf.defDim(ncid,'time',348);

dimid(2) = netcdf.defDim(ncid,'lon',144);

dimid(3) = netcdf.defDim(ncid,'lat',91);

[Madeleine P.]

Walter Roberson, your experiment seems to have worked! I'm now outputting netcdf files with spi. Thanks so much to you and Geoff for your suggestions and persistence in helping me! (even when I was about to give up). Cheers, M.P.

Answer 2

[Madeleine P.]

Didn't come up with an answer for this - added in the time for loop, reran the script, takes 6 hours to run, and I still end up with an empty array, Z[] (!). No error messages.

?????

[Geoff]

I don't think you paid heed to what Walter said in the comments... I have edited my answer accordingly.