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?

1 Answer 1

0

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)]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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