Discussion: Karttapullautin error

Hi all
Really don't know where else to turn, so trying here.

I've got my hands on some Lidar data and am using it with Karttapullautin. It creates the contour files fine, but no matter how small a piece of data I use, it always throws an error on the Vegetation generation step:
"GD Warning: product of memory allocation multiplication would exceed INT_MAX, failing operation gracefully.
gdImageCreate error at GD/Image.pm line 86"
It then fails on any subsequent step that requires that file as input, and I don't get any .png files as output.

Anyone else had this issue. I really don't think this is an available memory issue, as I've even tried processing an extremely small tile (40mx40m) with virtually no data points.
I'm using the latest 64bit version of Karttapullautin, but the same error occurs on the 32bit version.

Pulling my hair our here. Any suggestions?

Thanks

Have you looked at the classifications with lasinfo or by opening the file in OL-Laser? (There are other ways, but those two are easy.)

For lasinfo, try a command like

lasinfo -i filename.laz -otxt
then look at the resulting txt file

In OL-Laser, open the file, wait for it to load, then look for an info button.

My guess is you have either only ground points (classification 2) or some other lack of points to use for the vegetation step. KP seems to work with type 2 (ground) and everything else an "other" classification, so I'm interested to see what the problem turns out to be.

{I can't tell where you are, so my first guess is that your data is already in meters...}

If your data is in feet, you should convert it to meters with a reprojection step in las2las, or (not as good) use the xyz scaling factors in the pullauta.ini file. The "correct" factor is typically 0.3048 or the reciprocal 1/(0.3048). I prefer to reproject because it preserves the georeferencing. Scaling the inputs works but creates more work later.

It's counterintuitive, but if you give feet data to KP, it assumes it's meter data (it's just numbers to the algorithms). Typically, it means KP allocates space that's about ten times bigger than necessary. (It's that reciprocal factor, which is about 3.28, squared, so nominally 10.76) Older versions would just stop working, but that might not be true anymore because of the way KP chops up the input and does it in little chunks.

Also look at the reported ranges of the xyz data. You should be able to make sense of the units. Does the range of z make sense given the relief and tree or building heights? is it possible z is in cm or some other fraction of a meter?

If you have xy in m and z in cm, that would be a good use of a z-scaling factor in KP.

Where’s the lidar from, Rich? Can you share a link? Then maybe someone can outline the right steps to use the source with KP.

Hi guys

Thanks for the suggestions. I'm in Cape Town, so I think the data is in meters already, as that is the standard in South Africa. Also, the data definitely has more than just ground points.
The scaling factor sounds like an interesting idea. How do I know whether my data is in meters or in cms? Below is the output of lasinfo:

thanks again
richard


lasinfo (181108) report for 'Clovelly.las'
reporting all LAS header entries:
file signature: 'LASF'
file source ID: 0
global_encoding: 1
project ID GUID data 1-4: 00000000-0000-0000-0000-000000000000
version major.minor: 1.2
system identifier: 'OTHER'
generating software: 'LP360 from QCoherent Software '
file creation day/year: 337/2018
header size: 227
offset to point data: 512
number var. length records: 0
point data format: 1
point data record length: 28
number of point records: 15996178
number of points by return: 14398070 1462656 127236 8216 0
scale factor x y z: 0.001 0.001 0.001
offset x y z: 0 -3000000 0
min x y z: -55423.918 -3778271.888 -5.508
max x y z: -52167.508 -3776297.298 1314.740
the header is followed by 285 user-defined bytes
reporting minimum and maximum for all LAS point record entries ...
X -55423918 -52167508
Y -778271888 -776297298
Z -5508 1314740
intensity 1 4087
return_number 1 4
number_of_returns 1 4
edge_of_flight_line 0 1
scan_direction_flag 0 1
classification 1 12
scan_angle_rank -23 32
user_data 0 0
point_source_ID 11005 12028
gps_time 176671193.502007 176908020.138680
WARNING: there is coordinate resolution fluff (x10) in XY
number of first returns: 14398070
number of intermediate returns: 135448
number of last returns: 14397992
number of single returns: 12935332
overview over number of returns of given pulse: 12935332 2670932 357047 32867 0 0 0
histogram of classification of points:
6507875 unclassified (1)
5266051 ground (2)
2466 noise (7)
6823 water (9)
167 road surface (11)
4212796 overlap (12)

I think it fails because all coordinate values are negative. Before that bug gets fixed I'd suggest as a workaround trying to translate points enough to make them positive, something like:

las2las -i Clovelly.las -auto_reoffset -translate_xyz 60000 4000000 0 -olaz -o Clovelly.laz

... and translating processed data back, if needed.

Do you have any metadata files for this data? That would tell you what projection and units it is in and would help a lot to fix it.

Using Acme Mapper set to UTM coordinates, I got this for a fort in Capetown:

34H 262217 6243083 (The six-digit number is easting and the seven-digit is northing. The easting is 500,000 at 21 degrees East for zone 34.)

Honestly, it's worth your time to figure out the projection and reproject it into a convenient projection such as UTM Zone 34 (southern hemisphere), in my opinion.

min x y z: -55423.918 -3778271.888 -5.508
max x y z: -52167.508 -3776297.298 1314.740

x difference: 3000 meters, sounds reasonable
y difference: 2000 meters, sounds reasonable
z difference: 1300 meters, holy cow

Unless z is in cm, then just 130 m

I have no expertise in KP at this time.

13 m

Ah yes, I can't do math today.

A google maps search turned up a suburb of Cape Town called Clovelly.

It has a sea level to max elevation of ~500m (1600 feet), so the elevation numbers might be in feet. 3000x2000m would be reasonable for a single lidar tile. It does cover the hill mass in question here, assuming this is the place. There are some 1300m relief terrains on the coast SE of town.

I've seen lidar sets with mixed units before (meters for x,y - feet for elevation). Been a while though.

Here's a list of some common map projections in SA. Some of these mention "Axes: westing, southing (Y,X)." The negative numbers might be the easting and northing versions of these.

Great. Thanks all for the suggestions and ideas.
I don't have any meta-data files for this dataset, but I guess I can drop the City an email and pose any questions to them. I got the dataset from them and they seemed knowledgeable enough.

I will try re-project the dataset into a positive number realm, to see if that solves it, and I'll also investigate whether the z-values are in feet. I would have thought the z-range to be around 200-300m, so perhaps 1300 feet is correct. Again, I will ask the City about this.

Give me the weekend to play around with the data a bit, and hopefully I can reply with some results on Monday.

Thanks again.

Search the dataset folder for "proj", maybe. Typically metadata is in txt or xml format (which you can open with a text editor).

Karttapullautin kaputten

Well, the problem is not the negative values in the x/y coordinates. They all seem fine, and produce a contour file in the correct position. The range of z values is what seems out of wack and what seems to cause the problem.

To answer cedarcreek's question about the projection:
"The South African Datum referencing system is referred to as the Hartebeesthoek 94 Datum and uses the World Geodetic System 1984"

Projection: Transverse_Mercator
False_Easting: 0.0
False_Northing: 0.0
Central_Meridian: 19.0
Scale_Factor: 1.0
Latitude_Of_Origin: 0.0
Linear Unit: Meter (1.0)
Geographic Coordinate System: GCS_WGS_1984
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_Hartebeesthoek_1994
Spheroid: WGS_1984
Semimajor Axis: 6378137.0
Semiminor Axis: 6356752.314245179
Inverse Flattening: 298.257223563


So the next step is to follow your advice and try figure out how to convert this dataset to UTM 34H.

las2las -i Clovelly.las -o Clovelly_UTM.las --t_srs ???? --scale 0.01 0.01 0.01

Oh, and:
"In South Africa, we use a "2D + 1D" system, where there is a separate datum for horizontal and vertical co-ordinates. The vertical datum is sometimes referred to the Land Leveling Datum and is based upon the measurement of mean sea level (MSL) at tide gauges at Cape Town, Port Elizabeth, East London and Durban at the beginning of the 1900's."

Sounds like they're trying to make things difficult for me ...

a contour file in the correct position

I am quite sure there is bug. Even if it makes contour file it does not mean there is no bugs in vegetation or cliff or final file rendering code making it crash.

Or one might say it's not a bug. The AI of Karttapullautin doing the interpretation is young, sort of teenager, and still a bit sensitive. She just can't handle all the negativity you throw at it. Please be kind and more positive :)

Kids these days.

Rich: based on that las2las command, it looks like you're using libLAS's las2las. Do you also have the LAStools versions of those tools (las2las & las2txt) available to KP? libLAS's versions won't work, and may be what's causing the error in the first place?

I think I do have the correct version of las2las -- the command above was me googling how to do the conversion. At the time I didn't realize it was referring to a different product.
So far I've been (unsuccessfully) trying to work on getting the .las file converted from WGS84 Hartebeesthoek 94 to UTM34H, as I figured that this would solve the problem in KP. No luck, as mostly this is beyond me.
Now I'm going to go back to re-projecting it to get rid of the negative values and seeing whether I can make any progress that way.

Hi all (- and specifically Jarkko)

I'm really hoping that someone can see something that I can't.
I've processed my Lidar file to:
- only keep every 16th point
- remove the noise
- reclassify the ground points
- reproject it so that there are no negative values
- take just the portion of the map that I need

but still no luck. KP still throws an error on the vegetation step.

I'm making the .las file available here:
https://www.dropbox.com/sh/p0hjob710bpmseg/AABMMkS...
in the hope that someone with a bit more insight might be able to see what I can do to get this to work.

(Lidar and KP promises to keep orienteering alive in Cape Town - with no real mappers and limited areas, we're struggling to keep the sport going - which is why I'm really hoping to be able to sort this issue out)

Thanks again

y values are still negative. Have you tried hitting the translate command I suggested?

Translate did the trick. map

I used the translate command as follows:

las2las -i Clovelly.las -o Clovelly_Reproject.las -translate_xyz 60000 800000 0

which, after trimming the file gave me the file I uploaded to dropbox with the values:

scale factor x y z: 0.01 0.01 0.01
offset x y z: 0 -3000000 0
min x y z: 5776.06 -2977971.74 3.50
max x y z: 7581.63 -2976747.37 423.48
the header is followed by 285 user-defined bytes
reporting minimum and maximum for all LAS point record entries ...
X 577606 758163
Y 2202826 2325263
Z 350 42348

Although the header still reports negative values, everything seems positive in the actual data. What am I doing wrong?


And then how did you get the output back to the correct location after producing the files, so that it lines up to the correct location?

(At least it shows that it's user error rather than a system error!)

800000 doesn't translate enough. -3778271.888 + 800000 = -2978271.888

4000000 will give you positive figures.

Woohoo! It works! Thanks very much.

Last question then: how do I go about translating the .png and .dxf back to their actual locations? Is there a way to do this? Or do I just live with the georeferencing of the map being wrong?

Plenty of ways to do it. But I think in Ocad you can just edit real world coordinate offsets (menu->options-> scales).