[gdal-dev] Rotate UTM-Grid

Julian Zeidler gdal at zeidlers.de
Mon Jan 29 02:07:56 PST 2018


Hi Even,

Thank you for a terrific answer, as always. Worked like a charm
For others who might face the same Problem here the working python code 
for arbitrary rotation values:

from osgeo import gdal
from osgeo import osr
import math

Pixsize=1  #in units of coordinate system normaly meters
xsize=425 #size in pixel
ysize=335
epsgcode=25833
rotateAngle=5 #angle in degrees
ulx=385984.495969
uly=5818972.084099
filename='out.tif'
###### End config#######

dxx=math.cos(math.radians(rotateAngle))*Psize
dxy=-math.sin(math.radians(rotateAngle))*Psize

dyx=dxy
dyy=-dxx

gt=[ulx,dxx, dxy, uly,dyx,dyy]


ds = gdal.GetDriverByName('GTiff').Create(filename, 
xsize,ysize,1,gdal.GDT_Int32)


ds.SetGeoTransform(gt)
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsgcode)
ds.SetProjection(srs.ExportToWkt())
#delete and clean
ds = None


Best
Julian Zeidler


Am 24.01.2018 um 22:33 schrieb Even Rouault:
>
> On mercredi 24 janvier 2018 15:06:04 CET Julian Zeidler wrote:
>
> > Dear all,
>
> >
>
> > I have the challenge that i have to allign a UTM-GRID with a pysical
>
> > experiment.
>
> >
>
> > The model is rotated by 5 degrees compared to the UTM-Grid (EPSG:25833),
>
> > so that N is at 355°
>
> >
>
> > I have to rasterize multiple Vectors to a 1m resolution rotated-raster,
>
> > so that they can be pysicaly built and the results of a modell run
>
> > directly compared to the measured physical model.
>
> >
>
> > Is there any way to define such a custom rotated Projection so that i
>
> > could use it with gdalwarp and other tools?
>
> This can be done with gdalwarp but requires a preliminary step when 
> you create
>
> an initial blank target raster with the right dimensions, and attach a 
> geotransform that has
>
> rotation terms and a CRS. Then you can use it as the output of 
> gdalwarp in update
>
> mode.
>
> For example with the GDAL Python API;
>
> from osgeo import gdal
>
> from osgeo import osr
>
> ds = gdal.GetDriverByName('GTiff').Create('out.tif', width, height, 
> band_count)
>
> ds.SetGeoTransform(gt) # where gt is an array of 6 values with the top 
> left pixel coordinate, resolution and rotation terms.
>
> #See "Affine GeoTransform" in http://gdal.org/gdal_datamodel.html
>
> srs = osr.SpatialReference()
>
> srs.ImportFromEPSG(25833)
>
> ds.SetProjection(srs)
>
> ds = None
>
> Even
>
> -- 
>
> Spatialys - Geospatial professional services
>
> http://www.spatialys.com
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20180129/b5124cb3/attachment.html>


More information about the gdal-dev mailing list