Geoserver - GetFeatureInfo of raster/wms layer multiple points/location or bounding box

I am using Geoserver 2.8.1 version. I have netcdf and grib files, and I am uploading those to Geoserver through Geoserver extensions. After that, I use those data to show layer feature info in the map. (According to zoom scale and image size)I put arrows in the map to reference wind and wave directions. I calculate, to which angle the arrow will be positioned, by getting netcdf info like this;

localhost:8080/geoserver/wms?Service=WMS&version=1.3&bbox=23.75,32.75,24.25,33.25&layers=it.geosolutions:u10&query_layers=it.geosolutions:u10&FEATURE_COUNT=10&REQUEST=GetFeatureInfo&INFO_FORMAT=application/json&WIDTH=20&HEIGHT=20&x=10&y=10&time=2015-06-16T18:00:00.000Z

I make above request for each point changing only bbox(Lets say 1000 times average). It is killing my application. Is there any way that I can get feature info/values for all the cells/points in one request?

I know I can't do it with WMS - getFeatureInfo now, but I am open to any suggestion, only limitation is I need to use java.

Note: I am not looking for values of one point in multiple layers, I am looking for values of multiple points. I can already get values for multiple layers.

Solution

I solved the problem after trying many ways, but couldn't find the time to post the answer here until now.

First, I was developing my own WPS process. However, during searching for code examples, I discovered Geoserver already had a WPS process for getting original Raster data. So, I decided to try that out.

First important note before reading steps; If you are using SOAP Ui or sth similar to test requesting, instead of trying it with writing code. You should know that getting response from the request takes milliseconds for me, however Soap Ui takes 10-15 minutes to parse the returning json data. So, don't worry about the process being slow, it is fast actually.

Activate WPS in Geoserver through this guide
Then go to Demos -> WPS Request Builder
Choose Process as gs:RasterAsPointCollection
input raster should be the one you want. Select it from the combobox there. In my case WindWaveModel:u10
TargetCRS as yours; in my case EPSG:4326
interpolation must be nearest as you are trying to get data for each cell
Process outputs you must select application/json from combobox there. If you select any other GeoServer throws exception.
Don't fill emisphere (I think stands for hemisphere) and scale.
Optional: If you want to discover more, you can press DescribeProcess on the screen.
For the request body, Geoserver provides Generate XML from process inputs/outputs button at the bottom of the page. You must click this button for Geoserver to create the request XML. First, have a look at the XML and try to understand what it is.
First, you will realize, what you filled in the page form is there in the XML as input. Then, you will see there is more than what filled in the form; like the bounding box of your raster layer. It is automatically extracted from layer information and gs:RasterAsPointCollection process provides you with bounding box input. So that you can get only some part of your raster layer, instead of all of it.
Second, you will see there is process chaining, as the first <wps:input> is output of a wcs:GetCoverage service. Here was somewhat challenging for me, because I have time dimensions in my netcdf layers. And gs:RasterAsPointCollection was using wcs:1.1.1 specification. If you look at the wcs:1.1.1 specification you will see no time dimension input. There were so many try&fails for me there, I will skip the detail, as I will provide you with the final result in the XML form with the addition of time dimension input.

<?xml version="1.0" encoding="UTF-8"?><wps:Execute version="1.0.0" service="WPS" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.opengis.net/wps/1.0.0" xmlns:wfs="http://www.opengis.net/wfs" xmlns:wps="http://www.opengis.net/wps/1.0.0" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:wcs="http://www.opengis.net/wcs/1.1.1" xmlns:xlink="http://www.w3.org/1999/xlink" xsi:schemaLocation="http://www.opengis.net/wps/1.0.0 http://schemas.opengis.net/wps/1.0.0/wpsAll.xsd">
  <ows:Identifier>gs:RasterAsPointCollection</ows:Identifier>
  <wps:DataInputs>
    <wps:Input>
      <ows:Identifier>data</ows:Identifier>
      <wps:Reference mimeType="image/tiff" xlink:href="http://geoserver/wcs" method="POST">
        <wps:Body>
          <wcs:GetCoverage service="WCS" version="1.1.1">
            <ows:Identifier>yourWorkspace:yourRasterLayer(e.g. it.geosolutions:u10)</ows:Identifier>
            <wcs:DomainSubset>
              <ows:BoundingBox crs="http://www.opengis.net/gml/srs/epsg.xml#4326">
                <ows:LowerCorner>minLon minLat(e.g. 20.0 30.0)</ows:LowerCorner>
                <ows:UpperCorner>maxLon maxLat(e.g. 55.0 48.0)</ows:UpperCorner>
              </ows:BoundingBox>
              <wcs:TemporalSubset>
                <gml:TimePosition>yourTimeDimesion(e.g.2017-03-15T06:00:00.000Z)</gml:TimePosition>
              </wcs:TemporalSubset>
            </wcs:DomainSubset>
            <wcs:Output format="image/tiff"/>
          </wcs:GetCoverage>
        </wps:Body>
      </wps:Reference>
    </wps:Input>
    <wps:Input>
      <ows:Identifier>targetCRS</ows:Identifier>
      <wps:Data>
        <wps:LiteralData>EPSG:4326</wps:LiteralData>
      </wps:Data>
    </wps:Input>
    <wps:Input>
      <ows:Identifier>interpolation</ows:Identifier>
      <wps:Data>
        <wps:LiteralData>nearest</wps:LiteralData>
      </wps:Data>
    </wps:Input>
  </wps:DataInputs>
  <wps:ResponseForm>
    <wps:RawDataOutput mimeType="application/json">
      <ows:Identifier>result</ows:Identifier>
    </wps:RawDataOutput>
  </wps:ResponseForm>
</wps:Execute>

HttpClient client = new DefaultHttpClient();
HttpPost post = new HttpPost("http://localhost:8090/geoserver/ows?Service=WPS");
HttpEntity entity = new ByteArrayEntity(xmlString.getBytes("UTF-8"));
post.setEntity(entity);
HttpResponse response = client.execute(post);
// parse the json data in the response with a streaming model as it will be big

Right now, I am querying both u10 and v10 raster layers (making 2 WPS calls) for getting eastward_wind and northward_wind values. Then I calculate euclidian sum for each cell and in the end create a color map out of result values. I also calculate wind arrow directions in a similar geometric equation.

I hope this helps other people too.

Please let me know, if I can make this process faster, or optimize it any other way.