How to use the nansat.Domain function in nansat

meanLat = sum(lats, 0.0) / 4.

# get string with WKT representation of the border polygon
print 'BorderPolygon:', d.get_border_polygon(), '\n'

# longitude and latitude grids representing the full data grid
longitude, latitude = d.get_geolocation_grids()
print 'longitude shape: :', longitude.shape, '\n'

# make KML file with image borders (to be opened in Googe Earth)
d.write_kml(kmlFileName=oFileName + '03_preview.kml')

# -- Create Domain with stereographic projection and resolution 1000m
srsString = "+proj=stere +lon_0=5 +lat_0=60"
extentString = '-te -1000000 -1000000 1000000 1000000 -tr 1000 1000'
d = Domain(srsString, extentString)
d.write_map(oFileName + '02_stereo_map.png')
print 'Stereo Domain:', d, '\n'


print '\n*** domain_test completed successfully. Output files are found here:' + oFileName

''' Domain class is a container for geographical reference of a raster

A Domain object describes all attributes of geographical
reference of a raster:
    width and height, pixel size,
    relation between pixel/line coordinates and geographical coordinates,
    type of data projection (e.g. geographical or stereographic)

'''

# Create Domain object. It describes the desired grid of reprojected image:
# projection, resolution, size, etc. In this case it is geographic projection;
# -10 - 30 E, 50 - 70 W; 2000 x 2000 pixels
d = Domain("+proj=latlong +datum=WGS84 +ellps=WGS84 +no_defs",
           "-te 25 70 35 72 -ts 2000 2000")
d = Domain(4326, "-te 25 70 35 72 -ts 2000 2000")
d.write_map(oFileName + '01_latlong_map.png')
print 'Latlong Domain:', d, '\n'

# get shape
print 'shape:', d.shape(), '\n'

# Generate two vectors with values of lat/lon for the border of domain
lonVec, latVec = d.get_border()
print 'lonVec :', lonVec, '\n'
print 'latVec :', latVec, '\n'

# Get upwards azimuth direction of domain.
bearing_center = d.upwards_azimuth_direction()
print 'bearing_center :', bearing_center, '\n'

# Create domain with stereographic projection

def get_displacement_km(n1, x1, y1, n2, x2, y2, ll2km='domain'):
    ''' Find displacement in kilometers using Domain'''
    lon1, lat1 = n1.transform_points(x1, y1)
    lon2, lat2 = n2.transform_points(x2, y2)

    if ll2km == 'domain':
        d = Domain('+proj=stere +lon_0=%f +lat_0=%f +no_defs' % (lon1.mean(),
                                                                 lat1.mean()),
                    '-te -100000 -100000 100000 100000 -tr 1000 1000')
        x1d, y1d = d.transform_points(lon1, lat1, 1)
        x2d, y2d = d.transform_points(lon2, lat2, 1)
        dx = x2d - x1d
        dy = y1d - y2d
    elif ll2km == 'equirec':
        # U,V Equirectangular
        dlong = (lon1 - lon2)*np.pi/180;
        dlat  = (lat1 - lat2)*np.pi/180;
        slat  = (lat1 + lat2)*np.pi/180;
        p1 = (dlong)*np.cos(0.5*slat)
        p2 = (dlat)
        dx = 6371.000 * p1
        dy = 6371.000 * p2

'''Advanced processing of MODIS images:
        retrieve chl, tsm, doc with boreali
        generate images
        '''

        pnDefaults = {
            'lmchl': [0, 5, False],
            'lmtsm': [0, 3, False],
            'lmdoc': [0, 2, False],
            'lmmse': [1e-8, 1e-5, True]}

        borMinMax = [[pnDefaults['lmchl'][0], pnDefaults['lmchl'][1]],
                     [pnDefaults['lmtsm'][0], pnDefaults['lmtsm'][1]],
                     [pnDefaults['lmdoc'][0], pnDefaults['lmdoc'][1]]]

        dtsDomain = Domain(opts['srs'], opts['ext'])

        fileName = self.get_metadata('name')
        oBaseFileName = self.get_metadata('name').strip('"').strip("'")
        ncName = opts['oDir'] + oBaseFileName + '.nc'
        print ncName
        prodFileNames = {}
        for pn in opts['prods']:
            prodFileNames[pn] = '%s/%s.%s.png' % (opts['oDir'], oBaseFileName, pn)

        if os.path.exists(ncName):
            print '%s already exist!' % ncName
        else:
            # good bits for NRT
            #self.add_mask(cloudBits=[1, 4, 5, 6, 9, 10, 13, 15, 20, 21, 23, 28, 29, 30])

            try:

n.resize()

# make KML file with image borders (to be opened in Googe Earth)
n.write_kml(kmlFileName=oFileName + '_preview.kml')

# make image with map of the file location
n.write_map(oFileName + '_map.png')


# Make image reprojected onto map of Northern Europe
# 1. Create Domain object. It describes the desired grid of reprojected image:
# projection, resolution, size, etc. In this case it is geographic projection;
# -10 - 30 E, 50 - 70 W; 2000 x 2000 pixels
# 2. Reproject the Nansat object
# 3. Make simple image
dLatlong = Domain("+proj=latlong +datum=WGS84 +ellps=WGS84 +no_defs", "-te 25 70 35 72 -ts 2000 2000")
dLatlong.write_map(oFileName + '_latlong_map.png')
print 'Latlong Domain:', dLatlong
n.reproject(dLatlong)
n.write_figure(oFileName + '_pro_latlon.png')

# Reprojected image into stereographic projection
# 1. Cancel previous reprojection
# 2. Get corners of the image
# 3. Create Domain with stereographic projection, corner coordinates and resolution 1000m
# 4. Reproject with cubic interpolation
# 5. Write image
n.reproject() # 1.
lons, lats = n.get_corners() # 2.
meanLon = sum(lons, 0.0) / 4.
meanLat = sum(lats, 0.0) / 4.
srsString = "+proj=stere +lon_0=%f +lat_0=%f +k=1 +ellps=WGS84 +datum=WGS84 +no_defs" % (meanLon, meanLat)

n = Nansat(iFileName)

# List bands and georeference of the object
print n

# Write picture with map of the file location
n.write_map(oFileName + 'map.png')

# Write indexed picture with data from the first band
n.write_figure(oFileName + '.png', clim='hist')

# Reproject input image onto map of Norwegian Coast
# 1. Create domain describing the desired map
# 2. Transform the original satellite image
# 3. Write the transfromed image into RGB picture
dLatlong = Domain("+proj=latlong +datum=WGS84 +ellps=WGS84 +no_defs",
                  "-te 27 70.2 31 71.5 -ts 500 500")

n.reproject(dLatlong)
# Export projected satelite image into NetCDF format
n.export(oFileName + '.nc')

# Collect values from interactively drawn transect
# 1. draw transect interactively
# 2. plot the values
values, lonlat, pixlinCoord = n.get_transect()
plt.plot(lonlat[0], values[0], '.-'); plt.show()


# run tests of other nansat components
import test_domain
import test_nansat

n2 = Nansat('stere.tif')
n.reproject(n2)
n.write_figure(fileName=oFileName + '_proj_1onto2.png', bands=[1,2,3], clim='hist')

# Reproject onto grid of another file (destination file is in swath projection)
n.reproject()
n2.reproject(n)
n2.write_figure(fileName=oFileName + '_proj_2onto1.png', bands=[1,2,3], clim='hist')

# Reproject onto grids of lat/lon
dFromGrids = Domain(lon=lonGrid, lat=latGrid)
n2.reproject(dFromGrids)
n2.write_figure(fileName=oFileName + '_proj_on_grid.png', bands=[1,2,3], clim='hist')

# reproject onto automatically generated domain
dstDomainAuto = Domain(srs="+proj=latlong +datum=WGS84 +ellps=WGS84 +no_defs", ds=n.raw.dataset)
n.reproject(dstDomainAuto)
n.write_figure(fileName=oFileName + '_proj_1auto.png', bands=[1,2,3], clim='hist')

# export all data into NetCDF format
n.export(oFileName + '_0.nc')

# export one band to GeoTIFF
n.export_band(oFileName + '_2.tif', bandID=2, driver='GTiff')

# create new object from given domain and array
# 1. Reproject the current object
# 2. Get array with data
# 2. Create new Nansat object from the given array and for given domain
n.reproject(dStereo)
array = n[1]
nStereo = Nansat(domain=dStereo, array=array, parameters={'name': 'band1'})