0

I have a question that I assumed would be straightforward, but I'm unable to properly solve it. I have points on a rotated lat/lon grid, and I want to get their coordinates in another CRS. Currently, it seems that the only way I can manage to do this is via a double transform: (1) perform the opposite rotation to get to a regular lat/lon format, (2) move the regular lat/lon coordinates into the new CRS. An example:

from pyproj import Transformer

origin_lon = 10
origin_lat = 47
# Transform COSMO rotated grid coordinates into unrotated coordinates
rotated_to_unrotated = Transformer.from_crs( 
    '+proj=longlat',
    '+proj=ob_tran +o_proj=longlat +lon_0=-180 +o_lon_p={} +o_lat_p={}'\
    .format(-180+origin_lon, 90-origin_lat), 
    always_xy=True)

# Transform unrotated coordinates into the desired CRS
unrotated_to_CRS = Transformer.from_crs( 
    '+proj=longlat',
    'epsg:3035',
    always_xy=True)

# Perform the two rotations in successive order
coordinates=[(0,0), (1,1)]
transformed_coordinates=[]
for pt in rotated_to_unrotated.itransform(coordinates):
    transformed_coordinates.append(unrotated_to_CRS.transform(*pt) )
print(transformed_coordinates)
# >> [(4320999.999999999, 2654053.479931807), (4432566.859410212, 2765225.490478056)]

The coordinate transformation is correct and really fast; but it's just not very neat to have to use two coordinate transformations rather than one. How can this be done in a single step?

Improve this question

1 Answer 1

Reset to default
1

Right, I fixed it. It turns out I had the projection 'backwards'. The following code shows how to do it in one step:

from pyproj import Transformer

origin_lon = 10
origin_lat = 47
# Transform COSMO rotated grid coordinates into unrotated coordinates
rotated_to_unrotated = Transformer.from_crs( 
    '+proj=longlat',
    '+proj=ob_tran +o_proj=longlat +lon_0=-180 +o_lon_p={} +o_lat_p={}'\
    .format(-180+origin_lon, 90-origin_lat), 
    always_xy=True)

# Transform unrotated coordinates into the desired CRS
unrotated_to_CRS = Transformer.from_crs( 
    '+proj=longlat',
    'epsg:3035',
    always_xy=True)

# Perform the two rotations in successive order
coordinates=[(0,0), (1,1)]
transformed_coordinates=[]
transformed_coordinates2=[]
for pt in rotated_to_unrotated.itransform(coordinates):
    transformed_coordinates.append(unrotated_to_CRS.transform(*pt) )
print(transformed_coordinates)
# >> [(4320999.999999999, 2654053.479931807), (4432566.859410212, 2765225.490478056)]

# Transform COSMO rotated grid into the desired CRS
rotated_to_CRS = Transformer.from_crs(
    '+proj=ob_tran +o_proj=longlat +lon_0={} +o_lon_p=-180 +o_lat_p={}'.format(-180+origin_lon, 90+origin_lat),
    'epsg:3035',
    always_xy=True)
for pt in rotated_to_CRS.itransform(coordinates):
    transformed_coordinates2.append(pt)
print(transformed_coordinates2)
# >> [(4320999.999999999, 2654053.479931807), (4432566.859410212, 2765225.490478056)]
Improve this answer
2
  • Should I rotate COSMO_REA6 files, or use the RLAT RLON in COSMO_REA6_CONST_withOUTsponge.nc described under opendata.dwd.de/climate_environment/REA/COSMO_REA6/… ? More generally, do you know of a cosmo_rea6_util.py (without cdo) ? Thanks
    – denis
    Mar 8 at 10:21
  • Hi @denis, I don't know what your question is exactly, but I'd generally recommend using cdo where applicable, e.g. cdo remap{con/dis/nn}, as it works very well and is super fast. The rotation I perform here is purely for plotting purposes, the data itself is not touched.
    – Erik
    Mar 9 at 9:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.