I have a table that looks like:
Name
X
X
X
Y
Y
Z
Z
I have a polygon feature class with the attributes:
Place
X
Y
Z
Essentially, I want to join that table to the polygon, but I want each row in the table to be kept.
When I attempt a join, it only joins to one x, y, z and not all.
I am using ArcMap.
I have a SRTEXT like this: GEOGCS["My Data",DATUM["Not specified (.....
I want to get PROJ4TEXT of this STREXT.
Can I get it from QGIS. Because I want to insert it in Postgresql SPATIAL_REF_SYS table. (or I need PROJ4TEXT data to use this projection)
Answer
You can write the text (which is called WKT) into a text file mycrs.prj , and run gdalsrsinfo on that. Running in the OSGE4W (or Linux/MAC) shell
gdalsrsinfo mycrs.prj >out.txt
should report the proj4 string and the WKT definition for postgis.
I just started to play with leflet/geojson a little. But my coordinates are not rendered properly and I have no clue what is going on.
My coordinates are: 52.23943, 4.97599. They work correct with the setView function.
var map = L.map('leaflet_map').setView([52.23943, 4.97599], 15);
But using a geojasonFeature they are, hmmm, 'projected', somewhere east of Somalia.
var geojsonFeature = {
"type": "Feature",
"properties": {
"name": "Coors Field",
"amenity": "Baseball Stadium",
"popupContent": "This is where the Rockies play!"
},
"geometry": {
"type": "Point",
"coordinates": [52.23943, 4.97599]
}
};
var myLayer = L.geoJson().addTo(map);
myLayer.addData(geojsonFeature).bindPopup("I am a gjson point.");
Anyone who knows what is happening here?
EDIT
Out of pure curiosity I changed the coordinates around [4.976143930893815,52.23925499011473] and the point appears at its correct location. A known bug!?
Answer
I wouldn't call it a bug, just a matter of confusing and contradictory standards.
When talking about geographic locations, we usually use Lat-long. This has been codified in the ISO 6709 standard.
When dealing with Cartesian coordinate geometry, we generally use X-Y. Many GIS systems, work with a Geographic Location as a special case of a 2 D coordinate point, where the X represents the longitude and Y represents the Latitude. This order of coordinates, is exactly opposite that of the regular Lat-long notion.
Coming to your problem:
The map.setView takes a l.LatLong as an input, where the first cordinate is a Latitude, and the second is Longitude.
Hence when you want 52.23N, 4.97E, you pass in [52.23943, 4.97599]
The GeoJSON standard says that for any point, the first parameter is the X Coordinate (i.e. Longitude) and second parameter is the Y coordinate (i.e. Latitude);
Hence when you want 52.23N, 4.97E in GeoJSON, you need to pass [4.97599, 52.23943]
For further reading, go through this Q&A
I am running mcc-lidar on a win7 64bits with 16gb of ram. Nevertheless when I try to process a las file with 1, 800 000 points using -s 0.5 and -t 0.3 the program tells me:
Error: insufficient memory; if applicable use a larger post spacing.
I am able to process the data using -s 1.5 -t 0.3, so I figure the problem is not with the number of points.
The program is installed at c:\program files (x86), which I understand is for 32 bits programs. Is there some specific instructions to install it as a 64 bit application in order to use my 16 gb of RAM?
Answer
Currently there is not a 64-bit compile although the source code is available so, you can compile it yourself using a 64-bit compiler.
I find your issue odd. 1.8 million points is not a very large problem and <1m post-spacing is common. Usually vendors tile lidar into < 20 million points and MCC has never demonstrated memory issues on 15-20m points on 32-bit systems with as small as a 0.5 post pacing (s=0.25). Perhaps you have sparse post spacing over a large area and there is insufficient RAM to create a 0.5m spline surface.
Is your data, in fact, 1m post spacing? The starting scale parameter should be 1/2 the post spacing. This defines the resolution of the spline surface in the first model iteration.
One immediate solution would be to create smaller tiles. Another option is to thin the data by querying for last returns and running the model on the resulting data. If you do this you will need to play with the parameters.
Am trying to use python pyproj to do the above conversion (convert/transform from Hong Kong GRID to Hong Kong UTM), but am getting different values from the expected result.
My expected conversion result is: "North 836148, East 822745" to "North 828672.7449, East 12488562.2742". See my codes below for both zones (N and S):-
# https://epsg.io/102141
inProj = Proj(init='epsg:2326') # PCS: Hong_Kong_1980_Grid or EPSG:102140
outProj = Proj("+proj=utm +zone=49 +ellps=intl +units=m +no_defs") # PCS: Hong_Kong_1980_UTM_Zone_49N or EPSG:102141
# outPro = Proj(init='epsg:102141') # This returns RuntimeError
# Expected result: 828672.7449, 12488562.2742
x1, y1 = 822745, 836148
# x1, y1 = 836148, 822745
x2, y2 = transform(inProj, outProj, x1, y1)
print (x2, y2)
Or
# https://epsg.io/32749
inProj = Proj(init='epsg:2326') # PCS: Hong_Kong_1980_Grid or epsg:102140
# outProj = Proj("+proj=utm +zone=49 +south +datum=WGS84 +units=m +no_defs") # PCS: Hong_Kong_1980_UTM_Zone_49S or EPSG:32749
outProj = Proj(init='epsg:32749')
# Expected result: 828672.7449, 12488562.2742
x1, y1 = 822745, 836148
# x1, y1 = 836148, 822745
x2, y2 = transform(inProj, outProj, x1, y1)
print (x2, y2)
Could it be that my expected result is incorrect or I am missing something? I want to certisfy the expected output by going from: "Hong_Kong_1980_Grid_System" to "Hong_Kong_1980_UTM_Zone_49N / 49S"?
North East North1 East1
836148 822745 828672.7449 12488562.27
831905 822495 827480.7854 12473200.37
838208 816918 810095.1039 12481647.01
838127 816314 805922.5972 12482194.51
837452 820146 815300.3898 12484301.71
833764 816330 808333.7182 12478896.42
841476 830964 834944.8892 12495195.37
The result needs to match above table from Hong Kong 1980 grid to WGS84/UTM 49S
I've written a few Python Toolboxes (which are new at ArcGIS 10.1), but am yet to decide whether/when I should write them rather than Python Script Tools in a standard toolbox.
I thought the Online Help might enlighten me when it prefaces some dot points with:
Once created, tools in a Python toolbox provide many advantages
However, the five advantages listed all seem to be over not being able to use Python to write tools, and none seem to specify an advantage of Python Toolboxes over Python Script Tools.
The two advantages that I can think of are:
Am I overlooking the compelling case that led Esri to provide the Python Toolbox capability and, if so, what is it?
Answer
The two are very, very close in functionality but not completely equivalent.
Common to both
TBX, holding all the implementation in one PYT) and distribute via email or shared network drivesUnique to TBX files:
Unique to Python Toolboxes:
PYT over TBX.)A while ago when I was working on my first dozen or so PYT toolboxes I got flustered at how much of a hassle it was to set up a PYT for the first time, so I developed a tool called tbx2pyt. It'll take a TBX toolbox and convert it to a PYT with a minimal loss of code. In fact, the PYT that powers it was first a TBX. This may be a good way to transition existing tools to the Python Toolbox format if you so desire. At the very least, it makes it possible to set up your tools' parameters using the UI before switching to code.
We are building map apis for javascript integrated with our own map data server(tiles). (Though we can use openlayers or something else, but we can not for some thrid-part reasons).
In a word, we are re-writing something like openlayers.
The following is what we have done:
When a map is created, it will have the container(dom element),center(LatLng), and zoom . Then we will orginaize the tiles inside the container according to the center and zoom.
When user drag the container. We will record the increasement of the mouse position. Then caculate the new center for the map, and then re-caculate the required tiles and orginaize them.
Like this:
onMouseMove(dx,dy){
var dinx=caculateDegreeInCurrentResolution(dx);
var diny=caculateDegreeInCurrentResolution(dy);
this.center=new LatLng(this.center.lat+diny,this.center.lnt+dinx);
this.refresh();
}
It works as expected, but we found that it is hard to drag the map. It is not smooth like openlayers or other map library during the drag operation.
So I wonder if we are going the wrong way?
Answer
Have you considered looking at the OpenLayers source code to see how they do it? I can't say it's something I would ever consider doing as I'd rather just use OpenLayers as is. It has a very open license so there shouldn't be any issues there. If it's a techincal problem then fair enough, although make sure you've taken enough time to be absolutely sure you can't just use OpenLayers as rewriting a whole map api sounds like a nightmare. Apologies if you have alerady thought of everything but I'd being doing a disservice if I didn't at least mention it!
Back to the question at hand, you could look at the OpenLayers source code to see how they do it. I wouldn't say you should just rip it off but certainly take inspiration and see if there is anything obvious you are missing. I don't know how the licensing works for OpenLayers (and maybe someone more knowlegeable could enlighten me) but I'm sure that if you did use parts of the OpenLayers source you would just have to appropriately attribute it to them in your source code. Even if you don't you just should as good practice. Anyway, here's a link to the OpenLayers 'DragPan' module on GitHub - OpenLayers DragPan - You may need to do a fair bit of digging from there to find exactly what you are after, or you can start looking from the top - OpenLayers on GitHub
Hope that helps you in some way!
Using QGIS v. 1.9 (1.7 does not seem to be available for Ubuntu Linux 12.04), do the following:
Create new project, specify coordinate system by clicking on button in bottom right hand corner, enable on the fly CRS transformation, set projection to WGS 84 / UTM Zone 18S. Pass mouse over project window and UTM coordinates appear in coordinates window in status bar at the bottom.
Add delimited text layer in decimal degrees. When asked for projection, specify same projection as in (1). Points in the layer are plotted but, mousing over the map, the coordinates in the coordinates window are now in decimal degrees (even though the projection is set to UTM).
Is there something I am missing here? I have tried doing this in so many ways, so many times, but the result is always the same. I seem to be setting the projection correctly but only decimal degrees show in the coordinates window.
With thanks, Martin
Is there such a thing as a line/vertices symbology? If not a workaround? In some situations, particularly when editing, I would like to see my lines as well as all the vertices similar to the effect I get when using the Edit Vertices tool but on all features in the feature class not just the selected feature.
The solution should be dynamic, without having to create a separate point layer but having one line features class loaded twice in the TOC is acceptable.
Programmatic solution using an example (addin, ArcObjects & VB.Net) is acceptable.
I have the ArcEditor license in case there is a trick to use with representations.
Answer
If you have ArcEditor or ArcInfo license level you can your Cartograpic Representations to show vertices.
You can use maker (= point symbol) on top of a line to show all vertices. See option on vertices in ArcGIS Help:
In fact, when Sinnott published the haversine formula, computational precision was limited. Nowadays, JavaScript (and most modern computers & languages) use IEEE 754 64-bit floating-point numbers, which provide 15 significant figures of precision. With this precision, the simple spherical law of cosines formula (
cos c = cos a cos b + sin a sin b cos C) gives well-conditioned results down to distances as small as around 1 metre. In view of this it is probably worth, in most situations, using either the simpler law of cosines or the more accurate ellipsoidal Vincenty formula in preference to haversine! (bearing in mind notes below on the limitations in accuracy of the spherical model).
Source: http://www.movable-type.co.uk/scripts/latlong.html
What is the reason why law of cosines is more preferable?
Note: The quoted text has been updated by its author as mentioned below.
Answer
The problem is indicated by the word "well-conditioned." It's an issue of computer arithmetic, not mathematics.
Here are the basic facts to consider:
One radian on the earth spans almost 10^7 meters.
The cosine function for arguments x near 0 is approximately equal to 1 - x^2/2.
Double-precision floating point has about 15 decimal digits of precision.
Points (2) and (3) imply that when x is around one meter, or 10^-7 radians (point 1), almost all precision is lost: 1 - (10^-7)^2 = 1 - 10^-14 is a calculation in which the first 14 of the 15 significant digits all cancel, leaving just one digit to represent the result. Flipping this around (which is what the inverse cosine, "acos", does) means that computing acos for angles that correspond to meter-length distances cannot be done with any meaningful accuracy. (In certain bad cases the loss of precision gives a value where acos is not even defined, so the code will break down and give no answer, a nonsense answer, or crash the machine.) Similar considerations suggest you should avoid using the inverse cosine if distances less than a few hundred meters are involved, depending on how much precision you're willing to lose.
The role played by acos in the naive law-of-cosines formula is to convert an angle to a distance. That role is played by atan2 in the haversine formula. The tangent of a small angle x is approximately equal to x itself. Consequently the inverse tangent of a number, being approximately that number, is computed essentially with no loss in precision. This is why the haversine formula, although mathematically equivalent to the law of cosines formula, is far superior for small distances (on the order of 1 meter or less).
Here is a comparison of the two formulas using 100 random point-pairs on the globe (using Mathematica's double-precision calculations).
You can see that for distances less than about 0.5 meters, the two formulas diverge. Above 0.5 meters they tend to agree. To show how closely they agree, the next plot shows the ratios of the law of cosines:haversine results for another 100 random point pairs, with their latitudes and longitudes randomly differing by up to 5 meters.
This shows that the law of cosines formula is good to 3-4 decimal places once the distance exceeds 5-10 meters. The number of decimal places of accuracy increases quadratically; thus at 50-100 meters (one order of magnitude) you get 5-6 dp accuracy (two orders of magnitude); at 500-1000 meters you get 7-8 dp, etc.
I'm writing a routine for creating the convex hull envelop around shapefiles of islands. I could read the data, extract polygons vertices into a geometry and then compute the convex hull.
Now I need to save it to a shapefile of polygons.
Is there a way to transform the geometry (named convexHull in my code) returned by the convexHull method into a polygon geometry; or am I forced to parse all vertices of convexHull and create a ring of a polygon?
def doConvexHull(infile, outfile):
inH = ogr.Open(infile, 0)
if inH is None:
usage("Could not open file {0}. Exit.".format(infile))
layer = inH.GetLayer()
# get all polygons
thisGeometry = ogr.Geometry(ogr.wkbPoint)
for index in xrange(layer.GetFeatureCount()):
feature = layer.GetFeature(index)
geometry = feature.GetGeometryRef()
ring = geometry.GetGeometryRef(0)
points = ring.GetPointCount()
for p in xrange(points):
lon, lat, z = ring.GetPoint(p)
thisGeometry.AddPoint(lon, lat)
convexHull = thisGeometry.ConvexHull()
drv = ogr.GetDriverByName( "ESRI Shapefile" )
ds = drv.CreateDataSource( outfile )
if ds is None:
usage("Could not create file {0}".format( outfile) )
lyrname = "convexHull_${0}".format( layer.GetName() )
lyr = ds.CreateLayer( lyrname, layer.GetSpatialRef(), ogr.wkbPolygon )
thisFeature = ogr.Feature( layer.GetLayerDefn() )
thisFeature.SetGeometry( convexHull )
lyr.CreateFeature( thisFeature ): error from geometry
On my map canvas I have various point objects spread through several layers. My plan is to 'animate' those objects by a loop from start to a predefined resolution.
By now I'm just testing the algorithm with the vertical positions of the existing objects and an exportCanvas function which saves the canvas as a single image for each step.
Because this method isn't very handy I want to have the output of the loop directly to the canvas - delayed by e.g. 0.5 seconds per step.
Above is an example of the merged jpg-files to a gif
The problem is that my loop doesn't work very well at all and freezes my QGIS for the whole process time (+delay).
I know that I need some kind of hyperthreading or callback from QGIS to my plugin but at this point I have no clue how to transfer this idea to useful code
Down below: the code so far
# Animation: Button clicked
def btnAnimationClicked(self):
rate = 24
# Werte zum Färben der xyz
for i in range(0,rate,1):
self.caller(i, self.callback)
# Animation: Caller
def caller(self, val, func):
func(val)
# Animation: Callback
def callback(self, i):
yhi = 5699049523 + 400000
ylo = 5696625163
delta = yhi - ylo
rate = 24
ste = delta / rate
ubound = ylo + (ste * (i+1))
lbound = ylo + (ste * i)
values = (
('hi', ubound + 1, yhi, 'red'),
('is', lbound, ubound, 'green'),
('lo', ylo, lbound - 1, 'yellow'),
)
# HaKasten abgestuft färben
for layer in QgsMapLayerRegistry.instance().mapLayers().values():
if layer.name()[0:9] == "xyz_______":
ranges = []
for label, lower, upper, color in values:
symbol = QgsSymbolV2.defaultSymbol(layer.geometryType())
symbol.setColor(QColor(color))
rng = QgsRendererRangeV2(lower, upper, symbol, label)
ranges.append(rng)
attribut = 'y'
renderer = QgsGraduatedSymbolRendererV2(attribut, ranges)
layer.setRendererV2(renderer)
for layer in qgis.utils.iface.mapCanvas().layers():
layer.triggerRepaint()
#time.sleep(0.5) # optional
self.exportCanvas()
I have a layer of street segments that are within a given distance of a fire station. Now I want to make a border around these street segments. If I use Buffer, it puts a buffer around each segment. If I use Convex Hull it makes a border around the street segments but does not follow the street segments but rather cuts across them. What I really need is an outside buffer of the street segments that follows the outermost street segments.
Is it possible in QGIS?
The image below shows the problem. The lines in green represent street segments that can be reached from a given point within 4 minutes. The yellow is 6 minutes and the blue is 8 minutes from the same point.
When I put a convex hull around the green segments (shown in purple), the hull includes some of the yellow segments because a convex hull is based on the smallest polygon that contains all of the lines. I see that a concave hull is what I really need, but QGIS does not support the concave hull.
I have been pulling my hair out over this issue and am now wondering if it is a bug rather than something I am doing wrong.
The Problem: I have a vector layer that users can edit. Each feature on the vector layer has an attribute called "Status" and I am using a unique values renderer to style the features. The layer can contain point, line, and polygon geometries. I have successfully set the layer styling and features are styled according to their status attribute. However, when a user selects a feature to modify the geometry the editing handles are not displayed which means the user is unable to modify the geometry, but the style clearly changes to the edit style. Here are the relevant parts of my code to create the layer:
// Capture the styles for the default render intents of a Vector layer
var defaultSymboliser = OpenLayers.Util.applyDefaults({}, OpenLayers.Feature.Vector.style["default"]);
var selectSymboliser = OpenLayers.Util.applyDefaults({ strokeColor: "red", fillColor: "blue" }, OpenLayers.Feature.Vector.style["select"]);
var temporarySymboliser = OpenLayers.Util.applyDefaults({ strokeColor: "red", fillColor: "blue" }, OpenLayers.Feature.Vector.style["temporary"]);
var editLayerStyleMap = new OpenLayers.StyleMap(
{
'default': defaultSymboliser,
'select': selectSymboliser,
'temporary': temporarySymboliser
}
);
// Create a symboliser lookup that will be used to stye features in their default intent
var symbolisersLookup = {
'Draft': { fillColor: '#0000ff', fillOpacity: 0.4, strokeColor: '#0000ff', pointRadius: 6 },
'Submitted': { fillColor: '#00ff00', fillOpacity: 0.4, strokeColor: '#00ff00', pointRadius: 6 },
'Review': { fillColor: '#ff0000', fillOpacity: 0.4, strokeColor: '#ff0000', pointRadius: 6 },
'Processed': { fillColor: '#9900ff', fillOpacity: 0.4, strokeColor: '#9900ff', pointRadius: 6 }
};
// Add the symboliser and lookup rules to the StyleMap object, using the default intent
editLayerStyleMap.addUniqueValueRules('default', 'Status', symbolisersLookup);
editLayerRefreshStrategy = new OpenLayers.Strategy.Refresh({force: true});
editLayer = new OpenLayers.Layer.Vector("User Submitted Information", {
strategies: [new OpenLayers.Strategy.BBOX(), editLayerRefreshStrategy],
protocol: new OpenLayers.Protocol.HTTP({
url: "Services/UserSubmittedInformationHandler.ashx",
format: new OpenLayers.Format.GeoJSON()
}),
styleMap: editLayerStyleMap
});
mainMap.addLayer(editLayer);
If I comment out the line:
editLayerStyleMap.addUniqueValueRules('default', 'Status', symbolisersLookup);
Then the features are no longer styled according to their status value but the editing handles do now appear. So, is this an OpenLayers bug or have I not done something important that allows the select render intent to work correctly when a UniqueValueRules is applied?
Answer
I've tried everything I can think of but the behaviour of OpenLayers remains the same. I'm convinced this is a bug so will be logging an issue in the bug tracker.
In the meantime I have managed to implement a work-around, that whilst not ideal, will work until there is a fix for the issue. In the end the solution seems obvious but it took me some time to appreciate that fact! So, put simply I remove the StyleMap in the beforefeaturemodified event handler and then apply it again in the afterfeaturemodified event handler. The following is a very rough example of what I do:
editLayer.events.register("beforefeaturemodified", editLayer, function (event) {
editLayer.styleMap.styles["default"].rules = [];
// In addition it seems that virtual vertices are also not respected, their display property is set to "none"
// so we need to override that here.
modifyControl.virtualStyle.display = 'visible';
});
editLayer.events.register("afterfeaturemodified", editLayer, function(event) {
// Put logic here to save edits, etc
// Reapply the style map
addEditLayerStyleMap();
});
function addEditLayerStyleMap() {
var defaultSymboliser = OpenLayers.Util.applyDefaults({}, OpenLayers.Feature.Vector.style["default"]);
var selectSymboliser = OpenLayers.Util.applyDefaults({ strokeColor: "red", fillColor: "blue" }, OpenLayers.Feature.Vector.style["select"]);
var temporarySymboliser = OpenLayers.Util.applyDefaults({ strokeColor: "red", fillColor: "blue" }, OpenLayers.Feature.Vector.style["temporary"]);
var editLayerStyleMap = new OpenLayers.StyleMap(
{
'default': defaultSymboliser,
'select': selectSymboliser,
'temporary': temporarySymboliser
},
{
extendDefault: false
}
);
var symbolisersLookup = {
'Draft': { fillColor: '#0000ff', fillOpacity: 0.4, strokeColor: '#0000ff', pointRadius: 6 },
'Submitted': { fillColor: '#00ff00', fillOpacity: 0.4, strokeColor: '#00ff00', pointRadius: 6 },
'Review': { fillColor: '#ff0000', fillOpacity: 0.4, strokeColor: '#ff0000', pointRadius: 6 },
'Processed': { fillColor: '#9900ff', fillOpacity: 0.4, strokeColor: '#9900ff', pointRadius: 6 }
};
editLayerStyleMap.addUniqueValueRules('default', 'Status', symbolisersLookup);
// Apply the style map to the edit layer
editLayer.styleMap = editLayerStyleMap;
}
I'm going to have a look at the OpenLayers code to see if I can track down the bug and fix it.
BEGIN EDIT:
the line
modifyControl.virtualStyle.display = 'visible';
only gave me the virtual vertices at begin/end of each line segment, but no vertices in the middle (for inserting new vertices). I disabled the stylemap to find out the default values of the virtualStyle property and replaced the line above through:
modifyControl.virtualStyle = {
cursor : "inherit",
fillColor : "#ee9900",
fillOpacity : 0.3,
fontColor : "#000000",
hoverFillColor : "white",
hoverFillOpacity : 0.8,
hoverPointRadius : 1,
hoverPointUnit : "%",
hoverStrokeColor : "red",
hoverStrokeOpacity : 1,
hoverStrokeWidth : 0.2,
labelAlign : "cm",
labelOutlineColor : "white",
labelOutlineWidth : 3,
pointRadius : 6,
pointerEvents : "visiblePainted",
strokeColor : "#ee9900",
strokeDashstyle : "solid",
strokeLinecap : "round",
strokeOpacity : 0.3,
strokeWidth : 1
}
Is it possible to generate map tiles from QGIS?
I guess QGIS Server must have a way to do it, but is there a way to just generate the tiles and save them to disk?
With the new symbology renderer and scale-dependent rendering, QGIS is a powerful tool for generating dynamic maps.
Since I'm using the new-symbology renderer, qgis-mapnik is not an option, and since I'm using scale-dependent rendering, just generating a large georaster and using gdal2tiles to tile it is also not an option.
Answer
It is possible to generate tiles using python console, you can read about it here. Keep in mind, that you might need to modify the script.
However I would encourage you to use mapnik for generating tiles, as above mentioned qgis script did not work well for me.
Styling map is easy with xml for mapnik, you will not have any trouble with it. Some modifications for quantumnik will let you generate tiles directly from qgis.
I have installed Mapnik and successfully tested it. Now I am facing the problem of creating my own customized tiles. I have googled it but everywhere it was telling about OSM tile generation.
My vision is simple: I have a huge map data of my target area. I want to create my own custom map in Mapnik/Tilemill. Then create tiles and upload them into server. A step by step instruction would be helpful to me since I am a newbie in Mapnik.
The main problem: How can I create tiles with my own data?
Answer
Like some suggested... I think your path from now should be:
generate_tiles.py script to create your tiles with Mapnik in the bounding box and zoom levels you specify (to specify bbox, zoomlevels, modify generate_tiles.py, last few lines I think).This is just a rough outline that should get you started, you can ask more specific questions and I'll clarify.
I have a raster with float number values for example: 2.544459104537964 and I want to round these up to integers. I have seen other posts on ArcGIS like this one but when I use Raster Calculator in QGIS, Int(myraster + 0.5) only makes the first value (which is 0) to 0,5 and the last one (which is 100) to 100,5. The other pixels stay float.
Answer
This is a bit of a hack but it works: The QGIS raster calculator does not support rounding as far as I know but you can use GDAL to perform float to int type conversions.
Raster -> Conversion -> Translate
change the datatype from float to int by using the -ot option
gdal_translate -of GTiff -ot Int32 E:/float.tif E:/int.tif
I need to be able to measure distances between points, however the distance need to be calculated in relations to elevation. The points are homes of storytellers from the 19th century and the places which are mentioned in their stories. The distance therefore has to be "walking distance". A path along a valley will likely be shorter than a path over a mountain even though the actual straight distance is shorter. Attached is a screenshot illustrating my thinking. In the image, paths A and C would therefore be calculated shorter than path B. 
The points are from a CSV file but I also have a raster layer with the elevation data.
Answer
Achieving this goal is somewhat a basic task in GIS, however the method in QGIS might not be trivial. Your best chance is to use GRASS's r.walk function, which creates an anisotropic cost surface (dem+slope+other factors).
First, you have to create a friction surface as an input to r.walk. In your case it can be a single-valued raster (1.0) matching the extent of your DEM. You can create it with r.mapcalculator with the formula: A*0+1 where A is your DEM.
Next, you have to select a set of starting points from your CSV. These are the points, the accumulated cost surface will be calculated from. You have to create an individual cost surface from every starting point. It might be smart to define the ending points associated with every starting point in this step (in individual layers off course). After, you can run r.walk with the created inputs. The starting points can be in a single layer, you can iterate through them with the green arrow in the dialogue box.
Now in an ideal case, you have the cost surfaces and the ending points for every cost surface. In theory you could find the least cost paths with r.drain, but for me, it ended up in an error (python couldn't import the QgisRaster library). If you run in the same issue, you can go with SAGA's "Least cost paths" algorithm. It will create a point and a line layer for every ending point with the cost surface (use the iteration button again). After you have all the lines, you can merge them into a single shapefile with SAGA's "Merge shapes layers" tool.
This method can be very slow with the increment of points, so if you have a lot of them, you might try to automatize the method with python. It will be still a lot of time to calculate (especially the cost surfaces), but you don't have to create tons of ending point layers manually.
I have a table of attributes in which a field contains combinations of IDs from another shapefile (ex. '[6, 4]' where 6 and 4 are each individual IDs in my shapefile). Here's an example for visualization:
(the one in the left would be the aforementioned attribute table and the one in the right the shapefile).
I was wondering if there is some way I could select the features in my shapefile using the combinations in the separated table. Mostly I'd like to know if it's possible to Select by Expression using "what's selected in another table" as input.
Answer
yes it is possible, but it is a bit a wild ralley of expressions in my brain;-). Also possible, that there is a more ellegant and easier solution with python...Let's say your layer A contains your ids which should be selected by the pairs of values (field "link" in brackets, [1,3]) in your layer B. The first step is to make virtual field, e.g. "sel" in layer B with the expression shown in the image below. This will dynamically reflect the selection state of your layer B, which we will need later in layer A:
The next step is to build the expression in layer A to select the features which are related to the selected features in B. There we will use the aggregate function which is very powerful:
array_contains(
string_to_array(
aggregate( 'B',
'concatenate',
replace(replace(replace("link", '[',''),']',''),' ',''),
filter:="sel"=1,
concatenator:=', ') ) ,
to_string( "id" ))
The aggregate function collects all items from the link field, which are selected (filter:="sel"=1). then we have to replace the brackets [1, 2] and spaces with the 3 replace statements. it would be easier to get rid of them before. at the end we put the string back into an array and proof if the id from layer A will be within the array...
My model consists of a for loop, where the original raster is hillshaded from different angles. I want to combine the resulting images to a single multiband image.
What I already tried was to put the loop results together with Collect Values and then merge them with Composite bands tool. I found that it is possible to use result of Collect values as an input of Composite bands, but Collect values returns just one value, regardless of how many iterations there are, so the final output raster's bands are all the same. This is definitely not what I want.
Now I do it without any loop, but having eighteen or even more same tools in one model is not very elegant. If it was possible to do this with a loop (especially in the model builder), I'd like to know how.
EDIT: after Composite Bands, I run Principal Components, and the result is quite fine. I can do other analyses too (in other models), but I don't want to drop the model with Principal Components. While maintaining it, I don't like that my model looks like this:
My question is: how to replace the column of hillshade tools and output with a loop? Is it possible within Model Builder?
I have some line points representing manholes. I also have line features representing roads. Both of these features are currently living in a PostGIS database. I'd like to automatically move these points to the closest part of the line, so all of the points are arranged in a straight line.
Is there an appropriate PostGIS function that does this? Or, a function in QGIS?
Edit: I've found the "Snap Geometries to Layer" function in QGIS 3.0 works great for this, but I don't want to create another layer with the straightened points, I want to update the geometry of the existing points. Any way around this?
Answer
To permanently update the point geometries on the database level, run
UPDATE AS pt
SET = (
SELECT ST_ClosestPoint(ln., pt.)
FROM AS ln
-- WHERE pt. = ln.
ORDER BY pt. <-> ln.
LIMIT 1
);
in either the DB Manager in QGIS or pgAdmin (or the client of your choice).
This will find the closest line to each point and projects the closest point on that line to the given point as the new geometry, effectively 'snapping' them to the closest line.
(Note: in your picture, the point at the lower left will be snapped to the closest line, so to the one it kind of sits on already; if that is undesireable, you can specify an attribute of the line table to match any of the point's to only snap to those lines, see -- WHERE ... above; uncomment that line if you need or delete it).
This will make use of the spatial index on both
After the update, I suggest to reindex the table and run VACUUM ANALYZE to update the table statistics, as part of standard maintenance.
In Modelbuilder I have set up a model, where in one step I do a "Calculate Field" like "round((!OID!-1- %MyField% /15)/452,1)+1" . For the calculation i use the value of %MyField%, which actually is an Value of a specific Field in the first Row of a Table. I tried to retrieve is with "Get Field Value" - this works. But how do I ensure, that I have the value retrievet at the moment when I want to "Calculate Field"? Of cource I make an "Precondition" connection, but if "MyField" is 0 (Zero), the precondition says False, so the Calculate Field is not executed. To explain more: I iterate through features (A) in feature class. In the iteration model i create an fishnet (new feature class (B)) inside the feature. Then i add a field to the new feature class (B).Now i need to calculate this field using the formula mentioned above. In the formula i need to use a value from a field of the feature (A). Is there a way? (I use Arcgis 10.2)
As you see on the attached image, I labelled contour lines on the map (in meters). Which function will allow the "halo" to take the background color of the map?
On the image, the halo is greenish which diminish the look of the map.
I use ArcMAP 10.
What do I download from http://www.gisinternals.com/stable.php for Win7-32 to cut 1:1M map sheets out of Australia.ecw file in GDA94/Geographic and save them as geotiff in WGS84/Geographic? I would like to automate the process to cut all of the 1:1M map sheets then reproject and convert them to geotiffs.
(I am thinking gdalinfo, gdaltranslate, gdalwarp but I am not clear about the ecw support and don't understand the gisinternals release numbers.)
Answer
I suggest to visit http://www.gisinternals.com/stable.php. Depending on the Microsoft VC runtime version (2005 to 2013) and win32 vs win64bit, you get a full package of GDAL. Following the information link on the site, you get to know which formats are supported.
If you don't need mapserver binares, just ignore them. The installer will not set up a full mapserver automatically.
Regarding ECW support, some versions (like http://www.gisinternals.com/packageinfo.php?file=release-1400-gdal-1-11-mapserver-6-4.zip) include SDK 3.3 support, which enables you to write ECW files up to a certain limit without the need for a license.
http://www.gisinternals.com/packageinfo.php?file=release-1700-gdal-1-11-mapserver-6-4.zip comes with support for SDK 5.1 (meaning read-only usage), while the MSVC13 compiled version has no ECW support at all.
So choose the package that fits your needs. All have the standard GDAL and OGR utilities.
How can I use gdal_translate for AAIGrid so that the raster values are not more than 4 decimal places.
Answer
It can be done by adding -co DECIMAL_PRECISION=4 in gdal_translate command.
I have a line vector layer with a field name of "buffer_distance". I want create a dissolved buffer layer based on it's field. I write below code:
from qgis.utils import iface
from qgis.analysis import QgsGeometryAnalyzer
mc = iface.mapCanvas()
layer = mc.currentLayer()
Distance_Field = layer.fieldNameIndex("buffer_distance")
QgsGeometryAnalyzer().buffer(layer, "/tmp/buffer.shp", False, True, Distance_Field)
when run this code in QGIS 2.18 python console, It creates a new shapefile with all table fields of input layer but there is not any feature. What is problem?
Answer
You have several issues in your code.
So, following code works as expected:
from qgis.utils import iface
from qgis.analysis import QgsGeometryAnalyzer
mc = iface.mapCanvas()
layer = mc.currentLayer()
feat = layer.getFeatures().next()
att = feat.attribute('buf_dist')
QgsGeometryAnalyzer().buffer(layer,
"/tmp/buffer.shp",
att,
False,
False,
-1)
After running it at Python Console of QGIS, I got desired buffer (retrieved from "/tmp/buffer.shp"):
as it can see at above image.
Editing Note:
In this case I would use QgsGeometry class because it also has a buffer method. Following code does the work:
mc = iface.mapCanvas()
layer = mc.currentLayer()
feats = [ feat for feat in layer.getFeatures() ]
buffers = [ feat.geometry().buffer(feat.attribute('buf_dist'), -1).exportToWkt()
for feat in feats ]
epsg = layer.crs().postgisSrid()
uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"
mem_layer = QgsVectorLayer(uri,
'buffers',
'memory')
prov = mem_layer.dataProvider()
feats = [ QgsFeature() for i in range(len(buffers)) ]
for i, feat in enumerate(feats):
feat.setAttributes([i])
feat.setGeometry(QgsGeometry.fromWkt(buffers[i]))
prov.addFeatures(feats)
QgsMapLayerRegistry.instance().addMapLayer(mem_layer)
I tried it out with line layer of following image:
and it works as expected.
My goal is to create an image collection that has 100% free cloud over a small region of interest, for example a lake.
This script filters landsat 8 images based on location and cloud cover:
var nocloudimages = landsat8.filterBounds(ROI)
.filter(ee.Filter.lt('CLOUD_COVER', value))
.sort('system:time_start', true)
However, as we know, the 'CLOUD_COVER' accounted for the whole landsat 8 image percentage of cloud, not a particular Region of Interest (ROI) cloud cover.
Is there a method to achieve that?
Is it possible to programmatically (in python) access the attribute info of the last added feature in QGIS 2.0 before edits are saved? I'm failing to find a way to get the ID and subsequently other attributes of my feature of interest (the one most recently added).
idx = layer.fieldNameIndex('myField')
myList = []
lastEdited = ID THAT I'D LIKE TO RETRIEVE
for feature in layer.getFeatures():
if feature.id() == lastEdited:
print feature.attributes()[idx]
From this I'd like to get 'myField' of my most recent addition.
UPDATE: Following lead of this thread, I was able to achieve my goals via the Python console, but am failing to get it to work for my plugin.
def updateTableWidget(self, fid):
idx = layer.fieldNameIndex('myField')
myList = []
for feature in layer.getFeatures():
if feature.id() == fid:
myList.append(feature.attributes()[idx])
Via the python console, I was able to connect to my layer like so:
layer.featureAdded.connect(updateTableWidget)
After adding a QgsMessageLog to my function, the attributes were correctly identified. I now have a beginner python issue, however. Calling the function in my plugin as shown below gives me TypeError: updateTableWidget() takes exactly 2 arguments (1 given):
self.updateTableWidget()
The function was called from feature form python file.
Answer
To get the last added feature in QGIS I used the event AddedFeat:
self.connect(self.iface,SIGNAL("currentLayerChanged(QgsMapLayer *)") ,self.ChangedTarget) # in setup
self.connect(layer,SIGNAL("featureAdded(int)"),self.AddedFeat) # from the ChangedTarget
def AddedFeat(self,id):
# added feature by id.
feat = QgsFeature()
self.prelayer.featureAtId(id,feat)
This event has to be connected to the active layer so the connect/disconnect occurs in the ChangedTarget(self,layer): event.
This returns the feature as it is created, it wont have attributes yet but you can extract them based on a button click or other event qgsVectorLayer::committedAttributeValuesChange looks promising.
We have a 1-meter LIDAR DEM from a city.
A small subset can be downloaded from this link:
This screenshot shows raw DEM with gray palette (darker belts are streets, and greyish and whitish rectangles are blocks): 
This correspond to a place in Santo Domingo city, which can be seen in this Google screenshot: 
On average, blocks are "uplifted" ca. 2 meters from streets, which is not true. We want to have a clean DEM to generate stream network and topographic wetness index (TWI). With the DEM supplied (we don't have the original bands from laser scanner), the hydrgraphic network seemed to follow a rectangular layout, and TWI resulted in blocks pattern. These pictures show the results:
This is the stream network result, generated with r.watershed in Grass GIS: 
And this is TWI result, generated with SAGA: 
We tried some procedures to solve this inaccuracy without success:
1) Denoising tool. We applied r.denoise tool in Grass GIS, but got some problems with the module installation. We ran it again with a shell in Windows, and received an insufficient memory message.
2) Filters. We ran different types of filters (low-pass, median, mean, etc.), with different windows sizes, and trying to put weight in the direction of the streets (Grass GIS, SAGA, QGIS).
3) Geostatistics. We generated points cloud strictly over streets (tried 1000 and 2000 points), generated a variogram model and then ran an ordinary krigging to fill blocks. Variogram modelling, and ordinary krigging, was done in R, using different packages. We got a linear variogram, so we wouldn't rely on kriging results.
4) Other tools. Installed ALDPAT tool, but couldn't get it to work because the program couldn't read the DEM.
In all cases, results in terms of the drainage network wasn't good, because we couldn't avoid the rectangular stream network; also, TWI still resulted in a blocks pattern.
In particular, with OK interpolated result, we got a point-like pattern DEM that affected the network result. However, the effect of the blocks pattern was diminished.
Also, we took a look at this question and answers...
Filtering out canopies and buildings from DSM to have a bare earth elevation
...which redirected us to Whitebox Geospatial Analysis Tools, but couldn't convert our DEM to LAS format. Also, we were not sure about the effectiveness of Bare-Earth DEM tool for us, because it's designed for removing semi-transparent objects, not blocks incorrectly "uplifted", which is our case.
We still want to generate a high quality DEM to make our hydrographic analysis, but don't know what else can we try.
Answer
If you are open to using alternative software to solve your problem, then I can suggest the Remove Off-Terrain Objects tool of the cross-platform open-source GIS Whitebox Geospatial Analysis Tools (of which I am lead developer). I realize that you said in your question that you could not convert your data to LAS format, but the tool takes a raster, not LAS file, as the input. I think you were perhaps mixing this tool up with the Bare-Earth DEM tool, which does take a LAS input. You can import your raster DEM into Whitebox as a GeoTIFF file or an ArcGIS floating point binary raster (.flt) or any number of other common raster formats.
Here is another example of it's ability to remove buildings from a raster DEM:
Importantly, the algorithm is not a filter and therefore will not impact the elevations of grid cells outside of 'off-terrain objects' (buildings). This is critical if you want to use your bare-Earth DEM to calculate the wetness index. (That said, I'd question the usefulness of estimating the wetness index of an urban or dense suburban location. Furthermore, there is no stream network through the dense urban landscape you show in your image...I'm certain most of the streams have been culverted.)
EDIT
Actually, looking at the sample dataset that you uploaded, I don't know if your DEM is suited to use with the Remove Off-terrain Objects tool. I had thought that the image that you uploaded was simply suffering from poor symbology, but I see now that it is in fact terraced (i.e. there are large flat areas steps within the DEM). Take a look at the following profile from your DEM:
The tool relies on accurate measures of slope (which I don't think are possible with your DEM...this should be flagged up if you want to calculate TWI as well) and the preservation of sharp slope boundaries between the ground surface and the buildings. But in your case the DEM is severely smoothed and these sharp edges aren't as apparent (your houses also appear to be about 2-3 m in height, which is a bit odd). Are you able to get your hands on the original unprocessed LiDAR DEM, or better yet, the point cloud data? For your application, I would seriously consider re-interpolating the DEM.
I would like to create a square buffer from a point feature but I do not understand the code that goes into it.
Similar questions have been asked on the forums.esri website but that was over 10 years ago, and it did not work when I tried the code.
How do I create a square buffer from a point feature?
Answer
Try these steps with ArcMap 10:
For a Python solution:
Using SearchCursor and InsertCursor to create square buffers
Is there a list or any other recommendation of certain GPS hardware (esp. GPS USB sticks) working with QGIS? I have been off this topic for far over a decade, and i remember that in former days not all GPS devices worked well with QGIS.
Additional information: Preferred OS at this time is Windows 10.
Answer
So because neither a list nor any recommendations seem to exist, here is my (one-point item at this moment) gps device recommendations list representing the empirical research effort of last weekend.
The device mentioned in my comment above works absolutely well and straightforward:
1) Windows 10 recognizes and installs the device automatically. No driver setup needed:
2) My first disappointment, that this thing does not work at all was unfounded, because ...
At first sight I expected the cut edge of the label to flash, but here is only black plastique. Well, ok, the LED is not the lightest I've ever seen...
Ahead to the supplied GPSinfo software, selected COM4: port and 4800 BAUD (Ann.: BAUD is Bit/sec.), sat aquisition took less than 20 seconds:
And the last small step into QGIS GPS info window:
And the device was about €40.
I am looking to create polygons based on a point count, and hopefully find a way to automate this too. My example is, I have 1,000 points in my city boundary and I want to create polygons over the area of 100 people in the city boundary. So my result would be 10 polygons each covering a 100 people and not overlapping. I am using ArcGIS Desktop 10.3 with an Advanced level license.
I couldn't find the answer online.
I would want to automate this too so I could do this multiple times as the population grows.
I don't want to go by distance I want to go by point count alone.
I have a PostGIS query that will return several million rows:
SELECT
t1.id AS id1,
t2.id AS id2,
ABS(t1.mean_h - t2.mean_h) AS h_diff,
ST_Distance(t1.the_geom, t2.the_geom) AS dist
FROM tas_ponds as t1, tas_ponds as t2
WHERE
(t1.gid > t2.gid) AND
ST_DWithin(t1.the_geom, t2.the_geom, 17000)
When run in psql, I get an out of memory for query result error.
Googling suggests that this an error within psql rather than postgres/PostGIS. Would amending the query into the form SELECT ... INTO x FROM ... fix the problem? Are there any other recommended approaches for dealing with very large datasets?
Answer
Some poking around does confirm this is a Postgres client problem, independent of spatial or server considerations: the client has a limited amount of memory to buffer the results before displaying them on the screen, which you're exceeding.
The recommended approach to handle this is to use a DECLARE / FETCH approach to access the data in smaller blocks than the total result set. You could also create a view with components of the query (e.g. distance) to cut down on the memory needed for the query operation itself.
I have a given table in postgreSQL with OSM data which is used for routing purposes in Java with PGRouting. There exist two columns which I don't understand. They belong to osm2po. They are labeled with "flags" and "clazz". I know, that flags describe the kind of road. It means, whether it is driveable by car, bike,... The problem now is, that I do not understand the meanings of the values in the column clazz. I just know that they also describe the road type. Can anybody explain, what the values like 12, 15, 31,... mean? Does a documentation exist where the meaning of the values is written down?
Answer
The file osm2po.config, which can be obtained from the downloads tab on the osm2po page, contains a table with four column, defined as:
1) concurrent order
2) class (1-127)
3) default speed in kmh
4) allowed transportation type (optional) - since v4.5.30
And here are some sample rows, which I think explain where the 12, 51, etc, you are seeing comes from (in the second column).
wtr.tag.highway.motorway = 1, 11, 120, car
wtr.tag.highway.motorway_link = 1, 12, 30, car
wtr.tag.highway.service = 1, 51, 5, car|bike
wtr.tag.highway.living_street = 1, 63, 7, car|bike|foot
It does say that the class can be in the range 1-127, and there are only 23 different values in that config table, so I trust that covers all the ones you are seeing?
Following comment from OP, here are the the official osm highway docs
I have a few digital images that are not georeferenced. E.g.:
Img01.tif
Img02.tif
...
Img99.tif
On the other hand, I have a text file, which contains georeferencing information for these images. E.g., I have a file c:\tfws\list.txt which contains the following information:
Img01.tfw 0.25 0 0 - 0.25 456258,125 4569852,125
Img02.tfw 0.25 0 0 - 0.25 456586,125 4570001,125
Img03.tfw 0.25 0 0 - 0.25 456952,125 4570300,125
...
Img99.tfw 0.25 0 0 - 0.25 458412,125 4575123,125
I need to generate a TFW file for each image in the text file. E.g. create a file named C:\tfws\img??.tfw containing the following information:
0.25
0
0
-0.25
456258,125
4569852,125
Does anyone know of a tool able to generate these TFW files automatically?
Answer
It isn't nice code, but for your specific circumstance, here is some python that should save typing them out by hand:
f = open('list.txt', 'r')
lines = f.readlines()
for line in lines:
lineparts = line.split(' ')
outfile = open(lineparts[0], 'w')
for i in range(1, 4):
outfile.write(lineparts[i] + '\n')
outfile.write(lineparts[4] + lineparts[5] + '\n')
for i in range(6, 8):
outfile.write(lineparts[i] + '\n')
(real code would sanity check, but I assume you can check the results)
When I request a radius buffer in QGIS (Dufor). my buffer comes out as a 20-sided polygon rather than a true "circle". This has lead to some missed blocks when clipping underlying vector layers for population density analyses (e.g., they would be within a circle, but this 20-sided polygon's straight sides preclude fringe blocks from being included in the radius).
is there a way I can go into the radius buffer command lines and either increase the number if sides to make the area more circular, or make it truly a "radius" buffer?
Answer
The "segments to approximate" option is what you are looking for. So the number you put into that filed will be the number of sides/quarter. So the default is 5, so you end up with a 20 sided polygon. Put in 25, and you end up with a 100 sided on, so it becomes smoother.
I want to reproject a tiff raster file from projection A (ETRS89/UTM Zone 32) to projection B (DHDN/Gauss-Kruger zone 2). The two projections have a different datum so I have to use a datum transformation (NTv2 BETa2007 in my case). I want to create a new raster file with projection B because I want to use it in another software. So reprojection on-the-fly is not what I need.
I have tried Raster > Transform (or similar - I do not know the English menu) but I can not see how to define a datum transformation here.
How can I do a permanent reprojection of a raster file in Qgis with a datum transformation?
Answer
Gdalwarp is the tool to reproject, you find it in Qgis under Raster->Projektionen->Transformieren or standalone in OSGEO4W.
Basic command is
gdalwarp -s_srs EPSG:25832 -t_srs "+proj=tmerc +lat_0=0 +lon_0=6 +x_0=2500000 +y_0=0 +k=1.000000 +ellps=bessel +units=m +nadgrids=./BETA2007.gsb +wktext" input.tif output.tif
BETA2007.gsb should be in the same folder, or use absolute path to it.
In Qgis, select the loaded tif as input, and its CRS, and paste the following in the target-CRS:
+proj=tmerc +lat_0=0 +lon_0=6 +x_0=2500000 +y_0=0 +k=1.000000 +ellps=bessel +units=m +nadgrids=D:\path\to\your\BETA2007.gsb +wktext
I've been looking through all of the GRASS commands in the QGIS 2.10 processing toolbox and it seems like one of these must be able to do this but I haven't found one yet. There is 'v.sample' which can sample a DEM at a point location. Is there a way to do something similar to sample a DEM when calculating the length of a line segment?
I have many point features in close distance to each other and when I label them, the labels overlap. Is there a way that QGIS automatically seperates the lapels from each other and resolves the overlap?
Answer
This is a common cartographic problem. There's no one-size-fits-all solution.
As MrXsquared pointed out, by default QGIS will rearrange labels to prevent them from overlapping. If you have too many labels for them all to display without overlapping, QGIS will eliminate some of the labels. You can override that setting by choosing the option to "Show all labels for this layer (including colliding labels" - but obviously this creates the situation you're trying to avoid.
Reduce the number of labels. As Erik mentioned, colliding labels means you probably have too much information on your map. Use symbology to convey the same information.
Make the labels smaller, by:
Increase the label spacing, allowing labels to be placed further away from the feature they label. This method works if you have multiple labels for the same feature. See here for details: Colliding labels for point features in QGIS
Manually rearrange the labels using the data-defined placement and the Label toolbar. See this section of the QGIS user manual for instructions.
On QGIS 2.4.0-Chugiak the Openlayers plugin is found under Web menu rather than the Plugins menu. How can you use the Python console to add a Openlayers Google Maps layer?
The following code give us the error OpenLayers plugin not loaded.
try:
olplugin = qgis.utils.plugins['openlayers']
ol_gphyslayertype = olplugin.olLayerTypeRegistry.getById(0)
olplugin.addLayer(ol_gphyslayertype)
except KeyError:
print 'OpenLayers plugin not loaded.'
Answer
Using the direct way as you tried in your example looks much cleaner, however, I did not find out either how to achieve this in QGIS 2.4.
So long, as a workaround, you can navigate through the menu interface manually to open the desired map from the OpenLayers plugin. In the following code, set the variables mapProvider and openLayersMap to your desired entries from the QGIS->Web->OpenLayers plugin menu. See inline comments for further explanation.
mapProvider = 'Google Maps' #also use e.g. 'OpenStreetMap', 'Bing Maps' etc. as given in the Web->OpenLayers plugin menu
openLayersMap = 'Google Physical' #also use e.g. 'Google Streets', 'OpenStreetMap', 'Bing Road' etc. as given in the Web->OpenLayers plugin menu
webMenu = qgis.utils.iface.webMenu() #get object of the Web menu of QGIS
for webMenuItem in webMenu.actions(): #open the Web menu of QGIS and loop through the list of web plugins
if 'OpenLayers plugin' in webMenuItem.text(): #look for OpenLayers plugin entry in the Web menu
openLayersMenu = webMenuItem #and open it
for openLayersMenuItem in openLayersMenu.menu().actions(): #open the OpenLayers plugin menu entry and loop through the list of map providers
if mapProvider in openLayersMenuItem.text(): #if the desired mapProvider is found
mapProviderMenu = openLayersMenuItem #open its menu entry
for map in mapProviderMenu.menu().actions(): #loop through the list of maps for the opened provider
if openLayersMap in map.text(): #if the desired map entry is found
map.trigger() #click the entry to load the map as a layer
I have merged a bunch of shapefiles in ArcCatalog and get the "the selected file failed to draw...." error. Every thing has merged alright in the db.
I'm new to GIS and could use some help.
Each of the shapefiles that were imported to ArcCatalog are an aggregate of the following file extensions: •.dbf •.prj •.shp •.shp (xml) •.shx
Answer
There are some definite known limitations regarding using shapefiles, but if you are looking to consistently use large datasets with too many features for a shapefile, I'd see if merging the shapefiles into a feature class in a File Geodatabase would work with your workflow as geodatabases are much more efficient at handling large datasets. Additionally, as mentioned by @Michael_Miles-Stimson (see comment to other answer) there other additional benefits to using a geodatabase rather than a shapefile. The main consideration before committing to a switch away from the shapefile is to make sure all intended users of the data will be using software that can utilize the new data format. Shapefile is a very widely recognized standard accepted by most GIS software, though the File Geodatabase is usable by a growing number of GIS softwares as well and is definitely usable by ArcGIS & QGIS systems.
The following link documents some of the known limitations of Shapefile size that I think may be the cause of your issue: http://support.esri.com/em/knowledgebase/techarticles/detail/37447
I have a map of small remnants of forest along creeks for an entire county in California. I need to divide the whole map into equal sections so I can print it so you can actually see the sections of forest. I'd like a map with a grid which I can look at, and then look at a zoomed in map of (for example) section A1.
I've tried the Grids and Graticules wizard, which gives me a grid for the map, but when I zoom, the grid remains the same. I want the grid sections to be georeferenced with the map, so when I zoom in, the grid zooms also.
I'm using ArcMap 10.1.
I want to run the majority filter in ArcGIS Desktop with a window size of 4 or 5 i.e. greater than the default 3.
Is it possible, or is there some other command (or tool) which might help me achieve that?
Answer
A Focal Majority function does a very poor job at establishing a MMU. I would recommend using a sieve approach. This will provide an exact defined MMU.
I believe that GDAL has a sieve model and it is also available in our Gradient Metrics ArcGIS Toolbox. It is an easy procedure to implement. The ArcGIS steps for using sieve to establish a minimal mapping unit of 10 cells are, more or less, as follows:
MinCells = 10
tmp1 = RegionGroup(InRaster, "EIGHT", "WITHIN", "ADD_LINK", "")
query = "VALUE > " + minCells
tmp2 = ExtractByAttributes(tmp1, query)
outraster = Nibble(InRaster, tmp2)
My client would like his shapefiles reprojected and I would like to get a script that can batch process all the files located in the same folder. I am new to Python and ArcPy but am very interested. We are using ArcMap10 as our platform software. Any help or advice would be greatly appreciated.
I have a high-resolution WorldView-3 image that I am trying to orthorectify using
gdalwarp -rpc -to "RPC_DEM=dem.tif" input.tif output.tif
However the output is not aligned with Bing or Google maps at all,
I have tried various DEM sources both in WSG84 as GeoID format, but both give bad results and also tried different target EPSGs.
When I manually use the RPC model to orthorectify small tiles of the image by projecting the corners using code from https://github.com/gfacciol/IS18 the tiles do match Bing and Google maps very well. However this is a slow process and results in edges of the tiles not matching exactly. But it shows that the RPC model does contain the correct information to orthorectify the image.
This is the code that I use to orthorectify a tile within the image that does seem to work (even though it seems like a stupid solution):
from pyproj import Proj, transform
from rasterio.windows import Window
import utils # https://github.com/gfacciol/IS18/blob/master/utils.py
outProj = Proj(init='epsg:32629')
inProj = Proj(init='epsg:4326')
x1, y1 = transform(inProj,outProj,aoi['coordinates'][0][0][0], aoi['coordinates'][0][0][1])
z = srtm4.srtm4(lon_center, lat_center) # determines ellipsoidal height at AOI
myrpc = utils.rpc_from_geotiff(image)
x1, y1 = myrpc.projection(lon_upperleft, lat_upperleft, z)
x2, y2 = myrpc.projection(lon_bottomright, lat_bottomright, z)
easting, northing = transform(inProj,outProj, lon_upperleft, lat_upperleft)
easting2, northing2 = transform(outProj,inProj,lon_bottomright, lat_bottomright)
w, h = x2 - x1, y2 - y1
window = Window(x1, y1, w, h)
metadata = ds.meta
metadata.update(transform=rasterio.transform.from_origin(easting, northing, 0.3, 0.3))
metadata.update(width=w, height=h)
with rasterio.open('test_ortho.tif', 'w', **metadata) as dst:
dst.write_band(1, ds.read(window=window)[0,:,:])
This is the output I get with --config CPL_DEBUG ON:
gdalwarp --config CPL_DEBUG ON -rpc -to RPC_DEM="/home/ubuntu/SRTM/srtm_WSG_ellipsoid.tif" input.tif output.tif
GDAL: GDALOpen(input.tif, this=0x56001a5b0880) succeeds as GTiff.
MDReaderDigitalGlobe: RPB Filename: input.RPB
GDAL: GDALOpen(/home/ubuntu/SRTM/srtm_WSG_ellipsoid.tif, this=0x56001a5b69e0) succeeds as GTiff.
OGRCT: PROJ >= 4.8.0 features enabled
OGRCT: Using locale-safe proj version
OGRCT: Source: +proj=longlat +datum=WGS84 +no_defs
OGRCT: Target: +proj=longlat +datum=WGS84 +no_defs
GDAL: GDAL_CACHEMAX = 3069 MB
Creating output file that is 18490P x 14370L.
GDAL: GDALDriver::Create(GTiff,output.tif,18490,14370,1,UInt16,0)
Processing input file input.tif.
WARP: Copying metadata from first source to destination dataset
GDAL: GDALDefaultOverviews::OverviewScan()
GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=0,0,4061x4083 Dst=0,0,4622x3592
0..GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=4050,0,4112x4030 Dst=4622,0,4623x3592
.10GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=8153,0,2066x3996 Dst=9245,0,2311x3592
..GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=10204,0,2108x4122 Dst=11556,0,2311x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=12263,0,2065x4113 Dst=13867,0,2311x3592
20GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=14322,0,2062x4006 Dst=16178,0,2312x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=0,4025,4052x4117 Dst=0,3592,4622x3593
..30GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=4030,3996,2069x4121 Dst=4622,3592,2311x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=6085,3993,2073x4106 Dst=6933,3592,2312x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=8143,3969,2062x4136 Dst=9245,3592,2311x3593
.40GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=10201,3968,2126x4157 Dst=11556,3592,2311x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=12269,3978,2109x4274 Dst=13867,3592,2311x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=14322,3978,2062x4280 Dst=16178,3592,2312x3593
.50GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=0,8125,1981x4248 Dst=0,7185,2311x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=1966,8105,2065x4124 Dst=2311,7185,2311x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=4024,8097,2065x4138 Dst=4622,7185,2311x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=6084,8093,2061x4148 Dst=6933,7185,2312x3592
60GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=8141,8091,2062x4140 Dst=9245,7185,2311x3592
..GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=10199,8095,2078x4159 Dst=11556,7185,2311x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=12267,8124,2138x4208 Dst=13867,7185,2311x3592
70GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=14344,8229,2040x4158 Dst=16178,7185,2312x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=0,12224,2054x4160 Dst=0,10777,2311x3593
..GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=1967,12223,2143x4161 Dst=2311,10777,2311x3593
80GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=4023,12224,2140x4160 Dst=4622,10777,2311x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=6083,12226,2104x4158 Dst=6933,10777,2312x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=8142,12225,2108x4159 Dst=9245,10777,2311x3593
.90GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=10200,12227,2133x4157 Dst=11556,10777,2311x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=12264,12253,2113x4131 Dst=13867,10777,2311x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=14324,12303,2060x4081 Dst=16178,10777,2312x3593
.100 - done.
GDAL: GDALClose(/home/ubuntu/SRTM/srtm_WSG_ellipsoid.tif, this=0x56001a5b69e0)
GDAL: GDALClose(output.tif, this=0x56001a5fa700)
GDAL: GDALClose(input.tif, this=0x56001a5b0880)
This is in GDAL 2.2.3.
In GDAL 3.1.0dev I get an additional debug output (the rest is basically the same):
RPC: Short-circuiting coordinate transformation from DEM SRS to WGS 84 due to apparent nop
GDAL: GDAL_CACHEMAX = 3069 MB
RPC: Using GDALRPCTransformWholeLineWithDEM
Does anyone have an idea what is happening here?
I'd like to take this attribute value: "MAP_PAR_ID": "0160011000", and populate a new field with "16-11".
This works: concat(substr("MAP_PAR_ID",0,3), '-', substr("MAP_PAR_ID",5,3)) but i need to detect leading zeros and not include them
I've tried numerous variations of CASE ELSE statements but can't even detect leading zeros. Feeble attempt below:
CASE
WHEN substr("MAP_PAR_ID", 0, 1) = '0' THEN substr("MAP_PAR_ID", 1, 2)
ELSE substr("MAP_PAR_ID", 0, 3)
WHEN substr("MAP_PAR_ID", 5, 1) = '0' THEN substr("MAP_PAR_ID", 6, 2)
ELSE substr("MAP_PAR_ID",5, 3)
END
Would the python console make things a little easier?
Answer
You can use a regex replace for this:
regexp_replace('0160011000', '0+([1-9]+)0+([1-9]+)0+', '\\1-\\2')
regexp_replace("MAP_PAR_ID", '0+([1-9]+)0+([1-9]+)0+', '\\1-\\2')
The regex is broken down like this:
The regex_replace function is defined like this:
regexp_replace(string,regex,after)
So in the after section we use \\1-\\2 to add the values from group 1 and 2 into our new string.
I recently purchased a Dell XPS-8300 (i7-2600|16 GB DDR3 RAM|64-bit Windows7) which I intend to use primarily for GIS programs. I bought a refurbished model because it was much cheaper and because I am able to choose a graphics card of my own.
I am currently in the process of purchasing said card, but I haven't really seen much information about what kinds of cards are optimal for ArcMap or QGIS rendering as opposed to Battlefield 3 or Skyrim.
I am currently looking at an XFX AMD Radeon HD6870 PCIE 2GB Dual Mini or an XFX AMD Radeon HD 6870 900M 1 GB DDR5 DUAL MINIDP HDMI DUAL DVI PCI-E Video Card
The 2 GB model is about $40 more than the 1 GB.
I have a 460 watt power supply and do not know if the extra ram affects power consumption. The 1 GB model, from all reviews, should work fine with my system, so I am hoping the 2 GB would work fine as well.
However, since I already have 16 GB of onboard ram and have read that GPU RAM isn't really important unless you go under the amount needed for a given program and settings, I am wondering if this extra GB would be worth it for my purposes:
For example if I have a parcel dataset with hundreds of thousands of polygons, would that extra GB of GDDR5 RAM significantly improve the time it takes to render the polygons when I scroll/zoom in on ARCmap/QGIS; or does it mostly worthless due to the fact that I already have 16 GB that can be used?
Also are there any specific advantages/disadvantages to using a ATI vs NVIDIA card for GIS display programs in general?
Answer
Any old video card will work for the 2D display functions. The video card's 3D capabilities only come into play when using specific 3D GIS features such as ArcScene or ArcGlobe in ArcGIS Desktop. If you aren't planning on doing 3D visualization then it does not matter one bit.
I would spend the extra money on an SSD instead.
