Deformation when using reproject

[Juanes Marín]

Hi everyone, I'm developing a software that uses datasets from various radar locations in Colombia in order to estimate precipitation, previously I used ArcPy tools to process all the data, but now I'm trying to use opensource tools such as wradlib.

I'm having an issue when georeferencing the dataset, it comes in antenna coords from the radar itself (range,azimuth) but I need them to be in lat,lon WGS84. It seems like a starightforward process, but the result I'm getting looks distorted compared with the manual process.

The data should look circular-ish since the range is the same in all angles, but for some reason when I try to reproject it results in a elliptical shape, looks like this:

I've done a manual process transforming each pair of range,azimuh to lat lon and then make a shapefile of points wich gets interpolated using IDW obtainting this result:

The manual result covers the same extent as the radar should, for some reason the Y axis is exagerated when using georef.reproject. I've tried many geographical coordinated systems and always get the same result.

I'm new to posting questions in this kind of forums, please let me know if I missed something or broke any rule.

Here's the code that I'm using (there are a lot of commented and redundant lines because of all the testing that I've tried)

import wradlib as wrl
import matplotlib.pyplot as pl
import matplotlib as mpl
import numpy as np
import os
import pandas as pd
import netCDF4 as nc

def wradtest():

file_route =r"D:\Downloads\pyart_test2\CEM200916235306.nc"

with nc.Dataset(file_route) as ds:
data = ds['reflectivity'][:360].data
ran = ds['range'][:].data
az = np.linspace(0, 360, 361)[0:-1]
lat = ds['latitude'][:].data
lon = ds['longitude'][:].data

proj_wgs84 = wrl.georef.epsg_to_osr(4326)
proj_utm18n = wrl.georef.epsg_to_osr(32618)

radar_location = (lat, lon, 0)
sitecoords = (lat,lon)
polargrid = np.meshgrid(ran, az)

coords = wrl.georef.spherical_to_proj(polargrid[0],polargrid[1],0,radar_location,proj_wgs84)

#coords, rad = wrl.georef.spherical_to_xyz(polargrid[0], polargrid[1], 0, radar_location)
#geo_coords = wrl.georef.reproject(coords, projection_source=rad, projection_target=proj_wgs84)

x = coords[..., 0]
y = coords[..., 1]

xgrid = np.linspace(x.min(), x.max(), 2000)
ygrid = np.linspace(y.min(), y.max(), 2000)
grid_xy = np.meshgrid(xgrid, ygrid)
grid_xy = np.vstack((grid_xy[0].ravel(), grid_xy[1].ravel())).transpose()

xy = np.concatenate([x.ravel()[:, None], y.ravel()[:, None]], axis=1)
gridded = wrl.comp.togrid(
xy,
grid_xy,
max(ds['range'][:].data),
np.array([x.mean(), y.mean()]),
data.ravel(),
wrl.ipol.Nearest,
)
gridded = np.ma.masked_invalid(gridded).reshape((len(xgrid), len(ygrid)))

gridded[np.where(gridded == -9999.0)] = np.nan


with nc.Dataset(r"D:\Downloads\pyart_test2\200916235306_ELLI.nc", mode='w') as out_nc:
out_nc.createDimension('latitude', len(xgrid))
out_nc.createDimension('longitude', len(ygrid))
out_nc.createDimension('time', 1)

lat = out_nc.createVariable('latitude', "f4", ('latitude',))
# lat.units = radar.latitude['units']
# lat.long_name = radar.latitude['long_name']

lon = out_nc.createVariable('longitude', "f4", ('longitude',))
# lon.units = radar.longitude['units']
# lon.long_name = radar.longitude['long_name']

# tim = out_nc.createVariable('time_coverage_start', str, ('time',))

ref = out_nc.createVariable('DBZH', "f4", ('latitude', 'longitude'))
# ref.units = radar.fields['reflectivity']['units']
# ref.long_name = radar.fields['reflectivity']['long_name']

lat[:] = xgrid
lon[:] = ygrid

# str_out = np.array([radar.time['units'][14:]], dtype='object')
# tim[:] = str_out

ref[:, :] = gridded



if __name__ == '__main__':
wradtest()

[Kai Mühlbauer]

Hi Juanes,

wradlib uses (lon, lat) convention.

Please try to give sitecoords in (lon/lat), that might already fix the problem.

Best,
Kai

--
Diese Nachricht wurde von meinem Android-Gerät mit K-9 Mail gesendet.

[Juanes Marín]

Hi Kai,

It worked! Thank you so much. Don't know how I missed it.

--
You received this message because you are subscribed to the Google Groups "wradlib-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wradlib-user...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/wradlib-users/8403DA28-3847-428E-BAA3-BED41A495757%40uni-bonn.de.

[Kai Mühlbauer]

Hi Juanes,

Glad it works, and welcome to wradlib! Or let me rephrase, welcome to openradar community. Would be nice to see you around at https://openradar.discourse.group/.

This is the new place for asking questions and discussing matters of weather radar. And the good thing, Py-ART, BALTRAD and so forth are there too.

Best,
Kai

[Juanes Marín]

Thank you so much Kai, I'll surely be there from now on!

You received this message because you are subscribed to a topic in the Google Groups "wradlib-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/wradlib-users/WCvs0l5v6_o/unsubscribe.
To unsubscribe from this group and all its topics, send an email to wradlib-user...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/wradlib-users/425FA41A-3ECC-437B-A722-E1FDC5A4D2BC%40uni-bonn.de.

--

Juan Esteban Marín Ribón

Ing Civil