Nansat reproject fails for an ncep model wind field used for a SAR image over Svalbard

[mortenwh]

This should be tested using openwind software. Use the file RS2_20130520_062731_0076_SCWA_HHHV_SGF_260697_5426_8472420.zip for testing.

[asumak]

Solution: if gcps does not work well, use geoTransform. Can you please check e8bfb34, Morten? If it is ok, can you merge issue30 repository with the develop?

[akorosov]

Check usage of this option: -wo 'SAMPLE_STEPS=100'

[asumak]

SAMPLE_STEPS does not work well for this problem.

[knutfrode]

I have had similar problem before, with AVHRR images. The problem was that the GCPs are given as lon-lat (as almost always), and interpolation at high latitudes gives inaccuracies (singularity over pole).

The solution in that case was to transform all GCPs to polarstereographic coordinates, with the new function domain.reproject_GCPs() Try to run this function on the SAR image before warping the NCEP wind onto it.

[mortenwh]

Well - we think the problem is here related to the domain of the ncep data (0:360 deg in longitude). The region of the SAR coverage is then partly on the beginning/left and partly on the end/right side of a the ncep area. The solution we have come to is to reorganize the ncep array so that it covers e.g. -90:270 degrees instead... Asuka will implement this as soon as some other issues are solved, and then we'll see.

On 26 September 2013 10:04, knutfrode notifications@github.com wrote:

I have had similar problem before, with AVHRR images. The problem was that the GCPs are given as lon-lat (as almost always), and interpolation at high latitudes gives inaccuracies (singularity over pole).

The solution in that case was to transform all GCPs to polarstereographic coordinates, with the new function domain.reproject_GCPs() Try to run this function on the SAR image before warping the NCEP wind onto it.

— Reply to this email directly or view it on GitHubhttps://github.com/nansencenter/nansat/issues/30#issuecomment-25150099 .

[knutfrode]

Ok, that could also be the problem. But it should then be a generic problem, and thus it should be better to find a generic solution for all sorts of data, instead of making a fix for one specific mapper (NCEP). But didn't Anton solve exactly this problem a long time ago, at the early stage of Nansat?

[mortenwh]

I don't think Anton solved it already - he has been involved in the discussion on this issue... Anyway, I totally agree on a generic solution - that is important...

On 26 September 2013 10:56, knutfrode notifications@github.com wrote:

Ok, that could also be the problem. But it should then be a generic problem, and thus it should be better to find a generic solution for all sorts of data, instead of making a fix for one specific mapper (NCEP). But didn't Anton solve exactly this problem a long time ago, at the early stage of Nansat?

— Reply to this email directly or view it on GitHubhttps://github.com/nansencenter/nansat/issues/30#issuecomment-25152646 .

[asumak]

Nansat.reproject() gives warning if dstDomain is divided into two by reproject(). In 3e2c000, new Nansat object is created in nansat.shift(). the new object is shifted not to divide dstDomain. Currently we needs two steps as: n3 = n2.shift(n1) n3.reproject(n1) Do anyone have better solution?

[asumak]

Implementation Plan for create_shifted_vrt():

1. Shift geoTarnsform 1-1. if 'shirtDegree' < 0, then shirtDegree = 360 + shirtDegree 1-2. Shift geoTarnsform to 'shirtDegree' 1-3. if covered area < -360, then shifted +360. If covered area > 360, then shift -360
2. Switch images 2-1. Compute how many pixel want to shift using GeoTransform 2-2. Divide a band into two: 1st band : from 0, the size RasterXSize 2nd band : from RasterXSize- shirtPix 2-3. the 1st band starts from shirtPix and the 2nd starts from 0

Algorithm 1-1. if 'shirtDegree' < 0, then shirtDegree = 360 + shirtDegree 1-2. geoTarnsform[0] += shirtDegree 1-3. if geoTarnsform[0] < -360, then geoTarnsform[0] += 360 if geoTarnsform[0] + geoTarnsform[1] * RaterXSize > 360, then geoTarnsform[0] -= 360

2-1. shiftPix = shirtDegree / geoTarnsform[1] 2-2 & 2-3. modify VRT as : SrcRect xOff="0" xSize="720" DstRect xOff="shiftPix" xSize="720"

  SrcRect xOff="720-shiftPix" xSize="720"
  DstRect xOff="0" xSize="720"

[akorosov]

Closed in deb61a9f2f5176436e9175ad6340a81f0bf9ca0c deb61a9f2f5176436e9175ad6340a81f0bf9ca0c 4562f651745f4bc597946567c58af7306ad73357

VRT.get_shifted_vrt() is added, it rolls grib file (or any image which spans from 0 to 360) by -180 degrees (westwards). Rolling is implemented as described in #38

[knutfrode]

Does this really work?

When reprojecting NCEP wind over Radarsat2-images - I get 0-values for RS2-images crossing 0-meredian, but it works fine otherwise.

[knutfrode]

The problem is seen to be only for pixelfunctions. But again, the workaround of Alexander from #43 also works in this case:

After reading a (wind/grib) file with: n0 = Nansat(f)

One have to perform these two lines before reprojecting: n = Nansat(n0.vrt.fileName) n.fileName = n0.fileName

Then reprojection afterwards works fine also for pixelfunctions.

[knutfrode]

After updating my Nansat repository, I get now the following error when reprojecting a global NCEP wind field (0-360 East) onto a SAR image crossing 0 meridian:

File "/home/knutfd/software/nansat/vrt.py", line 1503, in get_shifted_vrt
    contents = node0.insert(cloneNode.rawxml(), 'VRTRasterBand', i)
TypeError: insert() takes exactly 2 arguments (4 given)

from the code:

w = Nansat(fw)
s = Nansat(fs)
s.resize(pixelsize=500)
w.reproject(s)

Before updating the repository, reprojection did not crash, but gave wrong values/directions in vicinity of 0-meridian.

[asumak]

insert() in nansat_tools.py was changed to the older version in 13d984c5fe. was it accidentally happened?

[akorosov]

No, not accidentally. I fixed errors introduced before. Adding childNode to the function call does not solve the problem of several layers of subNodes. What if we have to insert into child/childNode? Do we add childchildNode parameter to the function call? No, this is not the correct way to tackle such problems.

I forgot to fix get_shifted_vrt() to match the new insert() and I'll fix that.

[akorosov]

Fixed in f7c6a2b Node.replaceNode function is added to replace a node specified by a tag and absolute namber by a given node. Can be used with Node.insert which is also changed: it now returns Node and not XML

NB! VRT.set_subsetMask() also uses Node.insert and Node.replaceNode but is _NOT TESTED!_