geodesign / django-raster

Django-raster allows you to create tiled map services (TMS) and raster map algebra end points for web maps. It is Python-based, and requires GeoDjango with a PostGIS backend.
BSD 3-Clause "New" or "Revised" License
96 stars 39 forks source link

how to upload raster file and convert it? #28

Closed zaza316614 closed 6 years ago

zaza316614 commented 6 years ago

Hi.

I have to upload to server raster file and convert it when saving in file system if necessary. This is algorithm.

import sys, os, uuid, subprocess
from osgeo import gdal, gdalconst, osr

SUPPORTED_GDT = (gdalconst.GDT_Byte, )
SUPPORTED_GDT_NAMES = ', '.join([
    gdal.GetDataTypeName(i) for i in SUPPORTED_GDT
])

# os.chdir(r'C:\Users\Галина Геннадьевна\YandexDisk\Бизнес\ООО ГеоПризма\Разработка\Примеры данных\Тестовые данные\epsg3857 (pseudo mercator)')
os.chdir(r'C:\Users\Sergey\YandexDisk\Бизнес\ООО ГеоПризма\Разработка\Примеры данных\Тестовые данные\epsg3857 (pseudo mercator)')

input_file = "raster.tif"

gdal_source_ds = gdal.Open(input_file, gdalconst.GA_ReadOnly)
if not gdal_source_ds:
    sys.exit('GDAL library failed to open file {0}.'.format(input_file))

if gdal_source_ds.RasterCount not in (1, 3, 4):
    sys.exit('Only RGB, RGBA and single-band rasters are supported.')

for band_index in range(1, gdal_source_ds.RasterCount + 1):
    band = gdal_source_ds.GetRasterBand(band_index)
    if band.DataType not in SUPPORTED_GDT:
        sys.exit("Band #%(band)d has type '%(type)s', however only following types are supported: %(all_types)s."
                 % dict(band=band_index, type=gdal.GetDataTypeName(band.DataType), all_types=SUPPORTED_GDT_NAMES))

raster_xsize = gdal_source_ds.RasterXSize
raster_ysize = gdal_source_ds.RasterYSize
band_count = gdal_source_ds.RasterCount

print("InputRasterXSize=%(InputRasterXSize)d  InputRasterYSize=%(InputRasterYSize)d  InputBandCount=%(InputBandCount)d"
    % dict(InputRasterXSize=raster_xsize, InputRasterYSize=raster_ysize, InputBandCount=band_count))

# get the projection (SRS)
gdal_source_ds_projection = gdal_source_ds.GetProjection()
# provides the origin coordinates and pixel sizes, along with rotation values if the image isn’t situated so the top faces north
gdal_source_ds_transform = gdal_source_ds.GetGeoTransform()

if not gdal_source_ds_projection or not gdal_source_ds_transform:
    sys.exit('Only raster files with projection info are supported')

gdal_source_ds_osr = osr.SpatialReference()
gdal_source_ds_osr.ImportFromWkt(gdal_source_ds_projection)
gdal_target_ds_osr = osr.SpatialReference()
gdal_target_ds_osr.ImportFromEPSG(3857)  # 3857 - target spatial reference

reproject = not gdal_source_ds_osr.IsSame(gdal_target_ds_osr)

uuid = str(uuid.uuid4().hex)

# Separate in two folder levels by first id characters
levels = (uuid[0:2], uuid[2:4])
# path = os.path.join(os.getcwd(), 'tiled rasters', *levels)
path = os.path.join(r'c:\\', 'tiled rasters', *levels)

# Create folders if needed
if not os.path.isdir(path):
    os.makedirs(path)

output_file = os.path.join(path, str(uuid) + '.tif')

if reproject:
    cmd = ['gdalwarp',
           '-of', 'GTiff',
           '-t_srs', 'EPSG:%d' % 3857]
    if gdal_source_ds.RasterCount == 3:
        cmd.append('-dstalpha')
else:
    cmd = ['gdal_translate', '-of', 'GTiff']

cmd.extend(('-co', 'TILED=YES',
            '-co', 'BIGTIFF=YES',
            input_file, output_file))

subprocess.check_call(cmd)

gdal_target_ds = gdal.Open(output_file, gdalconst.GA_Update)
if not gdal_target_ds:
    sys.exit('GDAL library failed to open file {0}.'.format(output_file))

# Build overviews/pyramid layers
gdal_target_ds.BuildOverviews('average', [2, 4, 8, 16, 32])

raster_xsize = gdal_target_ds.RasterXSize
raster_ysize = gdal_target_ds.RasterYSize
band_count = gdal_target_ds.RasterCount

print("OutputRasterXSize=%(OutputRasterXSize)d  OutputRasterYSize=%(OutputRasterYSize)d  OutputBandCount=%(OutputBandCount)d"
    % dict(OutputRasterXSize=raster_xsize, OutputRasterYSize=raster_ysize, OutputBandCount=band_count))

del gdal_source_ds
del gdal_target_ds

How can I use django-raster to implement this?

Thank you.

yellowcap commented 6 years ago

It looks like you are reprojecting the raster to EPSG 3857 and then creating overviews within the raster file. Django-raster has a different approach, the raster is projected to 3857 but then it is split into many single tiles, which are then stored in the PostGIS backend that needs to be configured.

So you can not directly do what you have posted here. It would be possible to do with the GDALRaster class from Django, but the script would be similar to what you posted (except that it uses Django's internal gdal c-bindings and python api.

If you want to produce a raster you can display, you can use django-raster as described above. You would need to configure django with a postgis backend, as described in the docs.