Error in Wind correction function (broadcast shapes)

[simnh]

Hey again,

I get the following error:

generate_wind_correction - Start: 14:11:41:652020

Traceback (most recent call last):
  File "code/runme.py", line 22, in <module>
    generate_wind_correction(paths, param)
  File "/home/admin/projects/pyGRETA/code/lib/correction_functions.py", line 92, in generate_wind_correction
    A_cf_on = ((turbine_height_on / 50) * turbine_height_on / A_gradient_height) ** A_hellmann / resizem(Sigma, m_high, n_high)
ValueError: operands could not be broadcast together with shapes (1004,1043) (1200,1200) 

the m_high is 1200, however my matrix of all the other data has a different shape. Do you have any idea how to debug?

[kais-siala]

I find these dimensions (1004, 1043) a bit weird. Normally, all the maps are generated to cover the extent of the scope plus some pixels on the margins to have the same coverage as the low-resolution MERRA-2 data. So (1200, 1200) is more plausible to me. What is the geographic scope? I would like to try to reproduce your error.

[simnh]

I use a shapefile of Jordan, selected a box of 28 34 40 40 for the merra data. Or what specific information do you need

[kais-siala]

Hi, I have just downloaded a map of Jordan from GADM. I run the code (for generating the input files). I obtain maps with these dimensions:

[kais-siala]

As of MERRA-2: The map cardinal points are: W = 34.9576377899999997 S = 29.1858787500000005 E = 39.3020858799999999 N = 33.3681716899999969 If you follow the formulae provided in the documentation (instructions for downloading MERRA-2:

\begin{align*}
      minLat &= \left\lfloor\dfrac{s+0.25}{0.5}\right\rfloor \cdot 0.5 - \epsilon  \\
      maxLat &= \left\lceil\dfrac{n-0.25}{0.5}\right\rceil \cdot 0.5 + \epsilon \\
      minLon &= \left\lfloor\dfrac{w+0.3125}{0.625}\right\rfloor \cdot 0.625 - \epsilon \\
      maxLon &= \left\lceil\dfrac{e-0.3125}{0.625}\right\rceil \cdot 0.625 + \epsilon
\end{align*}

You obtain (with epsilon=0,01): minLat = 28.99 maxLat = 33.51 minLon = 34.99 maxLon = 39.385

So you cannot select any box containing Jordan. Otherwise you might have problems with the dimensions.

[simnh]

Thanks again, I updated the weather data but there is still the land use raster that has another size (same as before).

I think the problem is located at the landuse data I use. I have a landuse raster only for the area of the shapefile of jordan. This seems to cause the problem, I still need to map my landuse data to 1200x1200 pixel (of course the right ones).

[simnh]

I used this code in generate_landuse() now to create and read the (hopefully) correct subset raster using a window (LU_global points to my new landuse raster):

ds = rasterio.open(paths["LU_global"])
window = rasterio.windows.from_bounds(
        Crd_all[3], Crd_all[2], Crd_all[1], Crd_all[0], ds.transform, 1200, 1200)

with rasterio.open(paths["LU_global"]) as src:
     w = src.read(1, window=window)

At least the bounds look good for the new landuse:

BoundingBox(left=34.375, bottom=28.75, right=39.375, top=33.75)

The Broadcast error is solved, if it delivers the correct resulting raster.

I think if this is correct it could also be used for the code to be suitable with any landuse map that is not the world, but just a box bigger than the region of interest. Of course the pixels have to be general and not fixed at 1200. Let me know what you think

Sorry for the trouble, I just worked with rasterio and raster starting now with pygreta.