CRS of time series download

[David Pullinger]

I'm getting different "south_north/west_east" coordinates when I download a .nc file for either a bbox or a timeseries download for the same area.

When I download a bbox - the "south_north/west_east" appears correctly in EPSG:3035. However when I download a timeseries for the same area the "south_north/west_east" gives completely different (negative values), I've examined the .nc file of the time series and it includes the correct (desired) XLAT and XLON equivalent to the requested areas.

For comparison I've downloaded a box of data (for both timeseries and the same coordinates for the bbox) using the following code:

        parameter_url = 'https://wps.neweuropeanwindatlas.eu/api/mesoscale-ts/v1/get-variables'
        data_url = 'https://wps.neweuropeanwindatlas.eu/api/mesoscale-ts/v1/get-data-bbox'
        params = {
            'southBoundLatitude' : 50.61,
            'northBoundLatitude' : 51.0,
            'westBoundLongitude' : -4.17,
            'eastBoundLongitude' : -3.8,
            'height' : [height],
            'variable' : ['WS', 'WD'],
            'dt_start' : '2018-06-30T00:00:00',
            'dt_stop' : '2018-12-31T23:30:00'
            }

The coordinates that come out of this appear starting with: [-216000, -1296000] whereas they should be [3194160, 3325480] if they were in EPSG:3035.

Are the time series downloads in a different EPSG from the other sources? If so what is it?

 

[Neil Davis]

The timeseries data is defined on a lambert conformal grid. The WKT of the grid can be found in the crs_wkt attribute of the crs variable of the downloaded NetCDF file. It doesn't have an EPSG as it is a custom projection defined specifcally for the NEWA mesoscale simulations.

[Bjarke Tobias Olsen]

Hi David,

The projection, given as a proj4 string, is: '+proj=lcc +lat_1=30.0 +lat_2=60.0 +lat_0=54.0 +lon_0=15.0 +a=6370000 +b=6370000 +no_defs'

[David Pullinger]

Thank you both - that's resolved it.

[Dimitri]

Hi !

I have a question related to NEWA's XLAT/XLON parameters, which I think is relative to this post:

When I select an area and download time series for it (eg: t2, U100, Dir100, etc), I see that the .nc file contains 7 'south_north' and 7 'west_east' indexes that correspond (?) to 49 XLAT and 49 XLON parameters.

As, David, I want to assign ie: each of the 7 e.g: U100 time-series to a valid EPSG system.

I understand that I have 7 pairs (south_north, west_east) and expected to use XLAT/XLON parameters that come to 'readable' EPSG. But why do I have 49 XLATs, XLONs ?

Peculiarly, when I select on NEWA just a point and download time-series, then I have 1 pair of (south_north, west_east) and 1 pair of (XLAT,XLON).

I will appreciate any help !

[David Pullinger]

Hi Dimitri,

Is the problem you're seeing due to the change in the coordinate reference that is being used between the (south_north, west_east) indices when compared to (XLAT, XLON)? I've come across these issues when trying to import the NEWA download inputs into other software which requires a standard CRS (converting NEWA downloads into .wrg formats for example). (south_north, west_east) falls on a nice neat grid, with in your situation, 7 columns and 7 rows, but when XLAT/XLON are plotted you don't get neat columns and rows due to the warping resulting in 49 unique pairs (one for each point in the grid). It means that for coordinate transformation then I've had to loop through per point using something similar to below.

def ts_coordinate_converter(lon, lat):
    
    crs_wkt = 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6370000,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],PARAMETER["Latitude of false origin",54,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",15,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",30,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",60,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'

    custom_transformer = pyproj.Transformer.from_crs(
        pyproj.CRS.from_string(crs_wkt),
        pyproj.CRS.from_epsg(3035),
        always_xy=True
    )

    x, y = custom_transformer.transform(lon, lat)
    return (x, y)

 

[Dimitri]

Thanks, David

[Neil Davis]

The issue is that the time-series data is a raster in the projected space, but not in the lat/lon space. Therefore, the XLAT/XLON have to be two-dimensional, as the values at each point depends on both the x and y dimensions.

The approach from David above will give you a non-raster dataset in the new projection. If you want to get a raster back in the new projection, you can use the Windkit tool we have developed, which will "warp" the data from one projection to another, keeping the grid spacing approximately the same. When warping, you will be performing an interpolation, nearest neighbor by default.

import windkit as wk
import xarray as xr

ds = xr.open_dataset("mesotimeseries-Area 2.nc")
ds

ds_warped = wk.spatial.warp(ds, 3035)
ds_warped

Edited February 2 by Neil Davis
Fix code block