texts:
QGIS 2.8.1
Does anyone have the answer to why a generated heatmap appears without coordinates and therefore is not on top of the vectorlayer that it was generated from? I checked the CRS setting but that should be okay.
Any suggestions?


texts: QGIS - Drape (set Z value from raster) not working

texts: I have a shapefile layer containing lines and a raster layer imported from a Geosoft GXF.  I want to assign Z values to the nodes in the lines based on the intersected cells in the raster.  So I used the Processing Toolbox function Drape (set Z value from raster).  It doesn't work.  The Z values of the nodes don't change after running this function.  Is this thing broken?  Or is there something I need to know about the required inputs?  I'm displaying the vector and raster layers - the vectors are completely within the bounds of the raster.  They have the same CRS.


texts:
I'm looking to do the following.  Create a point on a map at a given coordinate (x,y) and then display text at that point.  Right now I can create a shape file but nothing is displayed on the map.  Hopefully someone can tell me why

Here is my code:

ArcDisplay.IRgbColor color = new RgbColorClass()
{
    Red   = 255,
    Green = 255,
    Blue  = 255
};

IPoint point = new PointClass()
{
    X = 0.0,
    Y = 0.0,
    Z = 0.0,
    SpatialReference = myMap.Map.SpatialReference
};

IElement elem = null;
ITextElement txtElement = new TextElementClass();
element = textElement as IElement;

ITextSymbol textSymbol = new TextSymbolClass();
textSymbol.Color = color;
textSymbol.Size = 16;
textSymbol.Font = (stdole.IFontDisp)OLE.GetIFontDispFromFont(new Font("Arial" , 16 , FontStyle.Bold);

element.Geometry = point as IGeometry;

textSymbol.Text = "My Test Text";
txtElement.Symbol = textSymbol;

((ITextElement)element).Symbol = textSymbol;

IGrpahicsLayer graphicsLayer = new CompositeGraphicsLayerClass();
((ILayer)graphicsLayer).Name = "Text Layer";
IGraphicsContainer graphicsContainer = graphicsLayer as IGraphicsContainer;
graphicsContainer.AddElement(element , 0);

myMap.AddLayer(graphicsLayer as ILayer);
((ILayer)graphicsLayer).Draw(esriDrawPhase.esriDPGeography , null , null);

Any ideas?


texts: MapBasic - Updating blank rows in column with the value above

texts:
I am writing a MapBasic tool to correct some incorrect values in a set of tables. So far I have managed to pull out the values I need in a new column:

What I need is to update any blank rows in the OBSERVATION column with the value above it (the L*** text).

How would I go about doing this in MapBasic?


texts: Exporting to SVG in ArcMap. How to change or call svg clippath polygon attributes in CSS

texts:
Recently I exported a map document from the ESRI ArcMap 10.2.2 to SVG.

I would like to animate this map with css and javascript however the "clippath" polygons that arcmap creates makes it impossible to do so.  How would can I get the group "Puerto Rico" to highlight over mouseover.

This is a summary of the code:

<svg width="1008pt" height="612pt" viewBox="0 0 1008 612" enable-background="new 0 0 1008 612" version="1.1" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
<g id="Puerto_Rico_and_Virgin_Islands">
    <g id="state" class="test">
        <g id="Puerto_Rico">
            <clipPath id="SVG_CP_1">
                <path d= coordinates are put in here />
            </clipPath>
        </g>
    </g>
</g>

<style type="text/css">
.test{
fill: grey;
stroke:#fff;
stroke-width:0.75;
stroke-opacity:1;
fill: orange;
}
#state:hover {
fill: red !important;
</style>

The real code:

visit jsfiddle for the entire code http://jsfiddle.net/jwitcoski/arv4g21f/3/

Using the wikipedia svg map that does not use clippath I was able to get everything to work.
http://jsfiddle.net/jwitcoski/y3yhgjjy/9/


texts: Identify the direction of traveling road segment

texts: I'm using pgRouting to find the least congested path using time as the cost metric. I am using average speed between two nodes to calculate time taken to travel in a segment. In pgRouting there is a forward cost and reverse cost. I have created tables to each and every edge to maintain reverse and forward speeds. How can I recognize whether a vehicle is travelling in forward direction or the reverse direction in sever side? I'm using postgreSQL as my database


texts: Labeling just a few specific features within a point dataset

texts:
I have a dataset with 95,000 points/features in and all I want to be able to do is label just a few specific features!

I have worked out how to label individual ones using 'Rule-based labeling' and using the syntax "Contact_ID"  = 'MBE_00618' (so that only the MBE_00618 target is labeled), however how do you alter the command so that it can label several targets in one go?

Not having any programming background am I struggling to find a solution.

The only way I have found for labeling multiple specific features is to create a new rule for each one, but is there a quicker, neater way?


texts: Find lat/long of Point 2 a given distance from Point 1

texts:
I'm creating a program using the Bing Maps API where the user inputs a bounding box with lat/long, and the program extracts a number of static images at zoom level 15 to later be mosaicked together.

If each square image is about 4100 meters wide, and the first image starts with its corner at (47.602030, -122.331572), what are the coordinates of the second image 4100 m to the east?

Here is an image to illustrate what I'm doing. The user-entered bounding box is black, the 4100 m long images are blue.

I'm using the RgoogleMaps package in R to access the Bing Maps API, so here is the R code used to calculate the ground distance of each image.

library(geosphere)
# bounding box coordinates (lower left , upper right lat/long)
extent1 <- c(47.602030, -122.331572, 47.619620, -122.291746)
# size of pixel in meters for zoom level 15 at given latitude
groundres <- (cos(extent1[1] * pi/180) * 2 * pi * 6378137) / (256 * 2^15)
# image height/width in meters for zoom level 15
img.size <- 640*groundres
# ground distance of bounding box (top corners)
extent1.m <- distRhumb(extent1[c(2,3)], extent1[c(4,3)], r=6378137)
# optimal number of images
ceiling(extent1.m/img.size)
# find lat/long of second image?

Edit: In addition to maptastik's answer, I found that package geosphere
has a function destPointRhumb() and similar functions that calculate a second point given an initial point, bearing, and distance.


texts:
I'm having some troubles trying to split a Linestring with 3 points.

My line is a circular route from a Garage (GA) to a pick up (PU) point, then a drop off (DO) point, and return to garage (img1)

And I'm trying to split this line by this three points, in order to get 3 lines (Garage(GA)->PickUp (PU), PU-DropOff(DO), DO-GA) First would be the left one, then the small one in the South, and finally the right one.

I'm using postgis ST_LineSubstring and ST_LineLocatePoint, but the result is not the expected, as I'm getting wrong return line:

I'm getting the PU-DO route (purple) as expected, and DO-GA as well (green) but for the fisrt segment (GA-PU) I'm not getting the one in the left, but I'm getting a new one behind the other (in blue)

I'm using the queries:

SELECT 1,ST_LineSubstring(geom, ST_LineLocatePoint(geom, PU),
 ST_LineLocatePoint(geom, DO))
FROM table;
SELECT 1,ST_LineSubstring(geom, ST_LineLocatePoint(geom, DO),
 ST_LineLocatePoint(geom, GA))
FROM table;
SELECT 1,ST_LineSubstring(geom, ST_LineLocatePoint(geom, GA),
 ST_LineLocatePoint(geom, PU))
FROM table;

And the results are the three lines shown.

Is there any way to get the other segment (the left one) in the third query instead of the one behind the others?


texts: QGIS and python: Keep a layer always on top of the others

texts:
I am working on a python plugin for QGIS using a layer containing linestrings, which I would like to have on top at all times. I have added the layer like the following:

QgsMapLayerRegistry.instance().addMapLayer(layer, addToLegend=False)

I don't want the layer to appear in the Layers Panel.

A similar question has been asked here without an accepted and for me meaningful answer.

According to Nathan W I can use the canvas method setLayerSet to set the render order on the canvas layers.

In case the current project already contains several layers, I need to know the current rendering order to correctly place my layer on top without disturbing the current order.
In case the user adds a new layer, how do I prevent the new layer from being ontop of mine? Is this handled by listening for some kind of a layer reordered signal on the canvas?

Question: What is the correct approach to handle this issue?

Update with solution notes:
I have implemented the answer provided by dmh126. To follow up on my questions I have determined that, it looks like canvas.layers() provide the layers in the order last given to the canvas.setLayerSet(...). Thus this list of layers can be used to include my layer on top directly:

layers = [my_layer] + canvas.layers()
canvas.setLayerSet([QgsMapCanvasLayer(x) for x in layers if x.id() != my_layer.id()])
# I don't want my layer in the list twice if it was already there due to a previous update.

The canvas signal layersChanged can be used to determine a change in the layers (new layer, removed layer, visibility change...), but remember to disconnect the listener when calling `canvas.setLayerSet(...)' because it will fire the signal again, probably causing a never ending call to the method called by the listener:

QObject.disconnect(canvas, SIGNAL('layersChanged()'), my_method)
layers = [my_layer] + canvas.layers()
canvas.setLayerSet([QgsMapCanvasLayer(x) for x in layers if x.id() != my_layer.id()])
QObject.connect(canvas, SIGNAL('layersChanged()'), my_method)

If I need to hide my layer, I can just remove it from the layer set like above:

layers = canvas.layers()
canvas.setLayerSet([QgsMapCanvasLayer(x) for x in layers if x.id() != my_layer.id()])


texts: Automatically style raster using unique values in QGIS?

texts:
In ArcMap, it is possible to automatically calculate the unique values for a raster and then apply a different style to each unique value (see 1st image).

However, in QGIS, I have to manually add values for styling when using the "Singleband pseudocolor" style (see 2nd image). Is there a way of auto-populating the unique values like ArcMap does?

ArcMap:

QGIS (how to auto-populate the value list with unique values?):


texts: How to merge linestrings which are connected with each other?

texts:
I would like to join linestrings which are connected to each other. The preferred output is a longer linestring based on all joined single linestings!
The table that holds my single linestrings is called “traffic_roads_temp”. Each linestring has the attribute ref.
Ref describes in my table the affiliation to a certain motorway or road.

As you can see in the picture below the motorway A 2 and A 7 could be easily joined by using something like this:

select ref, st_linemerge(st_collect(way)) from traffic_roads_temp group by ref

One of the problem is, that the output seems to be a multilinestring and not a linestring.

The main problem is the motorway A 39. It consists out of two separated road segments. Therefor I would like to generate two linestrings of the A 39 --> A 39.1 and A 39.2 All other motorways could be stored in one linestring per motorway.

Maybe somebody of the experts in this forum could help me to solve this problem(s):

Generate long linestrings from single smaller linestring. (No multilinestring)
I case of the motorway A 39 which consists out of two spread road parts – I would like to generate something like a series: A 39.1 and A 39.2

Problem-Sketch:


texts:
I downloaded the data in LAS file format and converted the LAS file to "multipoint" and then used "point to raster" tool to make it as raster, for missing data I used the raster calculator. I got the DEM but the problem is in the range of values (i.e, the lowest elevation and the highest elevation) which I am getting in my DEM, which I know that it is way more than what it should be, but not sure what are the factors which affect the range of the elevation values. What can be the possible reasons for this and how can I correct it (by looking for the solution I found somewhere that I should use some filter but I have no idea how to use a filter in this case)?

I want to keep the most of the elevation values as it, just want to remove the few points which are very different from the majority of elevation points and is responsible for making the range of values much bigger.


texts: Feature Layer Selection MVVM Silverlight

texts:
In my current project I have multiple feature layers which comes on Map when User Click on it. I'm able to represent them by MVVM on Map.
Now when User click /select by square area on layer, available geometry on it should be highlighted. I have to show the details of it on info Window. This data I need it in my viewModel for Highlight Grid Selected Rows Visibility too.

I am able to select the geometries by
https://developers.arcgis.com/silverlight/sample-code/start.htm#FeatureLayerSelection

But I m not getting the selected Data in my ViewModel & I m unable to show the infowindow by MVVM and Highlight selected Roads.


texts:
I'm trying to serve an existent MBTiles file using MapProxy. My YAML file looks like this:

services:
  wms:
    srs: ['EPSG:3785', 'EPSG:4326', 'EPSG:23028']
    image_formats: ['image/jpeg', 'image/png']

layers:
  - name: my_mbtiles_layer
    title: My MBTiles layer
    sources: [my_mbtiles_cache]

cache:
  my_mbtiles_cache:
    sources: []
    grids: [my_grid]
    cache:
      type: mbtiles
      filename: /path/to/myfile.mbtiles

grids:
  my_grid:
    base: GLOBAL_GEODETIC
    res_factor: 'sqrt2'
    bbox_srs: 'EPSG:4326'
    bbox: [-6.54,38.35,-3,40.65]

So, I defined my own grid, based on GLOBAL_GEODETIC, and used it for my MBTiles cache, that is the only source for my layer.

And the queries look like these:

http://myserver/service?LAYERS=my_mbtiles_layer&FORMAT=image%2Fpng&SRS=EPSG%3A4326&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&BBOX=-6.54,38.35,-3,40.65&WIDTH=256&HEIGHT=256

All my tiles form a coverage inside the BBOX (in epsg:4326): [-6.54,38.35,-3,40.65]

But the problem is I get a blank image of 256x256 px, and a HTTP code 200. So, the response is ok, but nothing is rendered in it. The MBTiles coverage is ok. I can see it in a desktop viewer. Here, an screenshot:

http://dl.dropbox.com/u/6599273/capturas/mbtiles_coverage.png

So, where can be the problem?

Many thanks in advance


texts: I want to read WaterML 2.0 time series. I am totally new to WaterML but I read it rely on GML 3.2. It would be great to use GeoTools to do that, maybe extending the GML classes. Is there some docs or examples about extending GeoTools GML? Is this a bad idea?


texts: Why is any pgr_* routing function taking forever based on OSM data in an pgrouting enabled DB

texts:
I loaded the German OSM dataset into the pgrouting DB by using osm2po 4.7.7. Everything works fine i have osm2po set up via it's config and it's working like a charm through it's Java part.

I had the *_2po_4pgr table imported without any problems. Even the *2po_v table gets imported, though I don't completely understand the relation of this table.

I executed the pgr_createTopology function which ran for quite a while (12000secs) while calculating all 6m edges. I thought this would do the deal, but still it's unbearably slow.

I would like to know if i forgot something. I was thinking of using pgRouting instead of the java library but at the moment its performance-wise just out of comparison.


texts: Adjusting z values for all vertices in polylines in ArcMap?

texts:
Is there way in ArcEditor (10.0) to adjust all the vertices z values by an offset? (translate in z axis?)

I have many polylines which need to be lowered (eg: subtract 1m from every vertex z-value).


texts: Avoid overwriting the style of a marked cell in Leaflet & Javascript

texts:
I have a system where I mark cells with different colours if a geoposition has been obtained within a cell. My issue is that due to the continuous update of geopositions, cells that were marked with one colour can be marked with the other colour as soon as the other colour's condition is true (else if in code). Here is my code

if (distance < positionerrortolerance){
    gpsPoint.properties.scanned ="Yes"
    gpsPoint.properties.distance = distance.toFixed(2);
    gpsPoint.properties.positionprecision = radius;
    rectangles.eachLayer(function (layer) {
        if (layer.complementId == k) {
          layer.setStyle({color: "#adff2f"});
        };
    });

 } else if (distance >positionerrortolerance && distance <sradius ) {
    gpsPoint.properties.scanned ="No"
    gpsPoint.properties.distance = distance.toFixed(2);
    gpsPoint.properties.positionprecision = radius;
    rectangles.eachLayer(function (layer) {
        if (layer.complementId == k) {
          layer.setStyle({color: "#ff7800"});
        };
    });
};

My idea is to assign after layer.setStyle({color: "#adff2f"}); an attribute = 1 to that layer and then check for rectangles.attribute != 1 in the else if (distance >positionerrortolerance && distance <sradius )I tried some things but they did not work.

Can somebody help? I assume there is an error when I assign the attribute.


texts: How to add a layer to SelectFeature without losing existing selections?

texts:
I'm using an OpenLayers.Control.SelectFeature for selections on multiple layers. However, when I add a layer using setLayer(), my selections on the other layers are lost.

Anyone know how to work around this? I would like to keep my existing selections on other layers when I add a layer to the SelectFeature control.

Here is an example:
MY EXAMPLE

Update:

I'm aware that this is part of the API. But i'm looking for a work around.

/**
 * APIMethod: setLayer
 * Attach a new layer to the control, overriding any existing layers.
 *
 * Parameters:
 * layers - Array of {<OpenLayers.Layer.Vector>} or a single
 *     {<OpenLayers.Layer.Vector>}
 */
setLayer: function(layers) {
    var isActive = this.active;
    this.unselectAll();
    this.deactivate();
    if(this.layers) {
        this.layer.destroy();
        this.layers = null;
    }
    this.initLayer(layers);
    this.handlers.feature.layer = this.layer;
    if (isActive) {
        this.activate();
    }
},


texts:
I am learning QGIS from scratch. Downloaded the software and data. I am following their manual uploaded in PDF format. I am stuck in the first step to add the vector layer. The layers panel show a black cross besides the layer name.

I believe it did not load the layer(S) properly. This is what I see on my screen.


texts:
I'm having some issues with QGIS 2.14 Essen (all of the versions...).
Specifically, I'm trying to do a buffer on a vector layer using the "Buffer Vectors" tool from the processing toolbox (ogr2ogr).

My layer is called "buffertest.shp", and contains 2 polygons. I'm using the default options, buffer distance of 1000 (Using CRS 27700, so it's 1000m), saving output as a temporary file. The console call is:

ogr2ogr.exe "\"[temporary file]\"" C:/Users/Skipper/Desktop/buffertest.shp buffertest -dialect sqlite -sql "SELECT ST_Buffer( geometry , 1000 ),* FROM 'buffertest' "

So far so good.

When the command is run, get a lovely message in the tool log saying:

cmd.exe /C ogr2ogr.exe "C:\Users\Skipper\AppData\Local\Temp\processing0c9922510f974f0eb4b34b3af7a0f9f6\9a2f675c807343c19e3a58fd66848382\OUTPUTLAYER.shp" C:\Users\Skipper\AppData\Local\Temp\processing0c9922510f974f0eb4b34b3af7a0f9f6\1476169842.153.shp 1476169842.153 -dialect sqlite -sql "SELECT ST_Buffer( geometry , 1000 ),* FROM 'buffertest' "
GDAL command output:
Warning 1: layer names ignored in combination with -sql.
ERROR 1: In ExecuteSQL(): sqlite3_prepare(SELECT ST_Buffer( geometry , 1000 ),* FROM 'buffertest' ):
no such table: buffertest
Converting outputs
Loading resulting layers

The following layers were not correctly generated.

Buffer

You can check the log messages to find more information about the execution of the algorithm

Which isn't ideal. No layer created, just the error there.
The log in the Log Messages Panel is the same, with the addition of QGIS saying it can't load the layer, because it doesn't exist.

Running the command through OSGEO4W shell works fine - gives the error:

Warning 1: layer names ignored in combination with -sql

But generates the right output.


texts:
I can't understand how to display all features from the feature layer. In the code bellow I specified 2 features in the layer and set in feature layer's mode to esri.layers.FeatureLayer.MODE_SNAPSHOT. According esri documentation this options means that all features will be selected and displayed. But they are not. The code is:

   var jsonFS = {
                "displayFieldName": "Name",
                "fieldAliases": {
                    "Name": "Name"
                },
                "geometryType": "esriGeometryPoint",
                "spatialReference": {
                    "wkid": 102100 //WGS_1984_Web_Mercator_Auxiliary_Sphere
                },
                "fields": [{
                    "name": "Name",
                    "type": "esriFieldTypeString",
                    "alias": "Name",
                    "length": 255
                }],
                "features": [{
                    "attributes": {
                        "Name": "bla"
                    },
                    "geometry": {
                        "x": -8919439.31450887,
                        "y": 4928270.761925456
                    }
                },
                {
                    "attributes": {
                        "Name": "bla bla"
                    },
                    "geometry": {
                        "x": -8155495.379532158,
                        "y": 5075380.311392084
                    }
                }]
            };
            var fs = new esri.tasks.FeatureSet(jsonFS);

            var featureCollection = {
                layerDefinition: {
                  "geometryType": "esriGeometryPoint",
                  "fields": [
                  {
                    "name": "Name",
                    "type": "esriFieldTypeString",
                    "alias": "Name"
                  }
                  ]
                },
                featureSet: fs
            };

        var featureLayer = new esri.layers.FeatureLayer(featureCollection, {
                mode: esri.layers.FeatureLayer.MODE_SNAPSHOT
            });


texts: How do I go about to perform a band ratio and Principal Component analysis in QGIS? I need to derive principal components from ratios of landsat TM bands .


texts: ArcGIS: How to create polyline features from the southern delineations of a polygon?

texts:
This can be considered a follow-up question to what I've asked a few weeks ago in this thread: Creating buffer only in specific direction using ArcGIS for Desktop?

Back then I was trying to set a buffer to several polygons in a specific orientation. The script david_p suggested worked perfectly, giving me results as shown in this picture:

What I need to do now is to measure the length of the buildings' southern-oriented walls (i.e. the shared lines of buffers and polygons). For this I want to create several polylines that should look something like the red line here:

I've already tried to create the polylines using the script from the link mentioned above with some adjustments. Unfortunately my knowledge of arcpy is very limited so that didn't work. Does anybody here know a way how I can create the polylines from the polygons without adding each single one manually? Due to the number of polygons and the fact that I'd have to do it for two different orientations this would take several hours of work.


texts:
I'm not too familiar with the GeoJSON functionality in PostGIS.

In essence - I'm having problems with the computation of ST_AsGeoJSON(ST_Union(g)) taking too long for large g.

Is it possible to instead store the corresponding GeoJSONs alongside their geometries in a table, and then directly merge the GeoJSONs to get the desired end product?

Would that be faster than first computing the union of geometries then converting it to GeoJSON?


texts: Get the intersection of two dataframes of coordinates R

texts:
I have two dataframes with geocodes. The first looks something like this:
spoints<- data.frame(x=c(1,2,3,4,5,6),y=c(6,5,4,3,2,1))

spoints maps a country.
My second data frame looks like this:
polyData<-data.frame(x=c(1,2,3,4,5,6,7,8,9,10),y=c(10,9,8,7,6,5,4,3,2,1),
col=c("a","b","c",etc.), id=c("a","b","c",etc.), average=c(44,33,66,55,etc))

This one contains coordinates to create voronoi clusters/polygons. but these are sketching polygons overlapping into the ocean.
So, I want to avoid that and have them stop at the country boundaries.
But now, I am having difficulty using GPC library or the others.
How can I do this?


texts: For a specific point in space, get the pixel in which that point lies (in R)


texts: With Python and ArcGIS, how do I add features to a featureclass?

texts: In ArcGIS Python code: CopyFeatures_management copies features to the same featureclass - that means overwrite to the feature class. How can I add new features to one featureclass or namely different featureclasses and prevent overwriting in the loop?


texts: Does QGIS use GPU’s to undertake calculations on Windows machine?

texts:
I am using SAGA Inverse Distance Weighted to interpolate between 2 large Raster files only it is akin many many hours to do 1 interpolation, let alone the 30 or so interpolations we have to do.
Does QGIS use GPU’s to undertake calculations on Windows machine or not as I am considering using a Cloud Service to try and speed up the process.
Alternatively, does anyone have any suggestions about how to speed up the process?  Please note, while it is interpolation, we are really extrapolating water surfaces.


texts: Is it possible to change the attribute that appears on the screen when I press the button "Identify features" above multiple points overlaid in QGIS?

texts:
In QGIS I need to map a lot of interest points in my city that was imported through a .csv file in order to compose the layer "Notas D6". I use the button "Identify features" to see more details about each point. But when there is some zone with overlaid points appears a list of numbers that isn't exist in attribute table of the layer. I don't know what are this numbers...

Attribute table of the layer:

Is it possible to change this "random numbers"(?) in order to select in this situations the  attribute that I prefer (corresponding with the column "Nota D6" in the attribute table shown above)?? This modification is important because allows a quickly identification of the point that I really want to consult more details instead of press in all numbers of the list to discover the point that I want.


texts: How to calculate the percentage of distribution in each polygon?

texts:
I have binary map for disease suitable area and I want to calculate the percentage using pixel count in each polygon how can I do that?

I use zonal histogram and do have frequency for each polygon is there is way to use this values to calculate the percentage in each polygon?


texts: Add points for every 100 m at polyline (QGIS/Mapinfo)

texts:
I have a polyline and want to set points for every 100 m. And if possible, with unique values for each point, 100, 200, 300 etc. Can it be done?

The QGIS-plugin, Sextante algoritmes "v.to points" and "v.segment" might just do that, but can figure it out.

Thanks


texts: Field Calculate with a function written in python within a larger python script

texts:
I'm writing a python program where I need to field calculate using a python function. The function takes two strings; one used for an if statement and the other to become an integer complete a simple addition equation. It then returns a string of that integer.

I've completed the calculation on the table in ArcCatalog and it works, but in Python it doesn't and gives me the following error:

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000989: Python syntax error: Parsing error SyntaxError: invalid syntax (line 1)
Failed to execute (CalculateField).

My function:

def yearStr(age, date):
    if age == 'BC':
         start = str((int(date) * (-1)) + 2251)
    if age == 'AD':
        start = str(int(date) + 2251)
    return start

Called using this:

yearStr(!Period!, !Year!)

In Python it looks like this:

arcpy.CalculateField_management(transp_table,"Start_Time","yearStr(!Period!, !Year!)","PYTHON_9.3","def yearStr(age, date):/n  if age == 'BC':/n    start = str((int(date) * (-1)) + 2251)/n  if age == 'AD':/n    start = str(int(date) + 2251)/n  return start")

I would appreciate any help! Have tried adding ' ' around field names and variables


texts: Hydrological triangle (geoscience) with QGIS 2.6.1

texts:
I want to interpolate ( TIN-Interpolation) between variable points (the best within a triangle) WITHOUT producing new shapefiles (or csv data file).

So what I have now is an interpolation of all points but I want to choose between which point QGIS is interpolating. Normally, I need to load a new layer with new (3) points (within the x,y,z- variable) everytime.

How can I choose between certain points that QGIS should interpolate?


texts: How to filter through geometry of feature collection instead of properties?

texts:
I have a FeatureCollection I made from a Polygon shapefile. I have fixed the geometries in qgis and should only be made out of many single polygons. But for some reasons some of my features are MultiPolygon and even some GeometryCollection which have linestrings and Polygons in the same features.
I would like to filter my features to only keep geometry type equal to Polygon. ee.Filter.eq only seems to filter through the properties of the features. is there a filtering tool to filter through other parts of the features?


texts: How to identify and highlight features programmatically?

texts:
In QGIS there is an option to select features of a selected layer using a "identify feature" tool. Then the feature get highlighted in a different color and attribute dialog is poping up.
like that i want to highlight the multiple features programatically.

Ex: I have a set of vector layers which consists of number of roads as features. I need to select few roads by looking at its attributes and highlight.

Is there any way of achieving this using python code?


texts:
New to this site... I'm working with Postcode regions in the UK and was wondering if anyone know if there is a dataset which give you postcode region and the most likely station destination within London...

i.e.  - LE - Leicester - St.Pancras / Kingscross


texts: Is MODIS sinusoidal tile boundary in Eathdata website real?

texts:
I was using MOIDS data MCD12Q1 (MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 500m SIN Grid V00). I try to use PyProj in Python to transform the Sinusoidal projection to lat/lon. It actually worked, but I mentioned that the the lat/lon boundary I have transformed is not consistent with the boundary information supplied in the Eathdata website.

Details:

tile h27v05, in file MCD12Q1.A2017001.h27v05.006.2018257173250.hdf, it own metadata:

UpperLeftPointMtrs=(10007554.677000,4447802.078667)
LowerRightMtrs=(11119505.196667,3335851.559000)

according to the mcd12 user guide, I set proj4 to
"+proj=sinu +a=6371007.181 +b=6371007.181 +units=m" then get the upperleft and lower right corners' lat/lon. But they are different from the boundary in the website. As shown below:

https://search.earthdata.nasa.gov/search/granules/granule-details?p=C186286578-LPDAAC_ECS&g=G1552253044-LPDAAC_ECS&m=22.11328125!108.703125!3!1!0!0%2C2&tl=1531459844!4!!&q=MCD12&ok=MCD12&sb=113.765625%2C23.765625%2C123.2578125%2C31.5703125&ac=true

Why do the tiles in the web site seem to be overlapping?


texts:
As the title suggests, how can I change the GPS tracker symbol in QGIS desktop?
Little explanation: I connected a GPS to my QGIS trough a COM port, the position is showed accurately on the map, but it has a default symbol as a marker.
I want to change the default (black circle with inner cross) to my own marker (a boat).

Haven't found it out yet how to do this.. Hope somebody knows this!


texts: pgShapeLoader: how to load a shapefile in a different schema than 'public'

texts:
I'm trying to load a shapefile in PostGIS using pgShapeLoader.

It works but I can load table only in the "public" schema and I don't know how to configure, the tool to upload my shapefile in another schema.

Is it possible?

Any suggestion for workaround?

Anyway I can use shp2pgsql and then psql from command line in PostGIS...


texts: hide attribute with null value using leaflet

texts:
http://jsbin.com/hinekaj/edit?html,output

this is a link to a map am working on. It will be great if someone could help me to hide the attributes that have "null" as their value. I am aware that I need to use the IF command but have been struggling with it.

        onEachFeature: function (feature, layer) {
          layer.bindPopup('Address: '+ feature.properties.address
                          + '

Sample 1: '+ feature.properties.sample_1 + ' Pb'
                          + "

Sample 1: "+ feature.properties.sample_1 + ' Pb'
                        + " , Sample 2: "+ feature.properties.sample_2 + ' Pb'
                         + "

Sample 3: "+ feature.properties.sample_3 + ' Pb'
                         + " , Sample 4: "+ feature.properties.sample_4 + ' Pb'
                         + '

Sample 5: '+ feature.properties.sample_5 + ' Pb'
                         + ' , Sample 6: '+ feature.properties.sample_6 + ' Pb'
                         + '

Sample 7: '+ feature.properties.sample_7 + ' Pb'
                         + ' , Sample 8: '+ feature.properties.sample_8 + ' Pb'
                         + '

Sample 9: '+ feature.properties.sample_9 + ' Pb'
                         + ' , Sample 10: '+ feature.properties.sample_10 + ' Pb'
                         + '

Sample 11: '+ feature.properties.sample_11 + ' Pb'
                         + ' , Sample 12: '+ feature.properties.sample_12 + ' Pb'
                         + '

Sample 13: '+ feature.properties.sample_13 + ' Pb'
                         + ' , Sample 14: '+ feature.properties.sample_14 + ' Pb'

         );
        },


texts: How to use label maximum line length for word wrapping in QGIS 3

texts:
I am confused on how to use the maximum line length for word wrapping a label in QGIS 3. The dialog and documentation show an option to "Wrap lines to" and then allow selecting "Maximum line length" but where do you specify the length?
I want to actually constrain the label to certain width (maybe even a rectangle) and have the lines split into multiple lines based upon the width.
Here is the documentation:


texts:
Trying to write a 3-band JPEG using gdal_merge but keep getting the following error:

Format driver JPEG does not support creation and piecewise writing.
Please select a format that does, such as GTiff (the default) or HFA (Erdas Imagine).

I'm not understanding this as on GDAL's Raster Formats page the Creation column for JPEG says 'Yes.' Additionally, the JPEG File Format help page indicates that it is possible to create JPEG.

Here is my command:

gdal_merge.py -of JPEG -co QUALITY=90 -co EXIF_THUMBNAIL=YES -n 0 -a_nodata 0 -o /path/to/Output.jpg /path/to/inputs/band1.tif /path/to/inputs/band2.tif /path/to/inputs/band3.tif

Why am I recieving this error? I realize that my input bands are 32-bit float but I just figured they would be truncated since JPEG only supports 8-bit. Could this be my problem? I don't want to spend time fixing that if there is something else I'm missing as to why I can't create jpeg.


texts: Render a map showing lots of individual traces/routes simultaneously

texts:
Google's "My places" is great for showing lots of traces (in my case, paths representing bike rides I've done). But it only shows 20 or so per page.

I'd really like a map (a statically rendered .png would be ok - but prefer dynamic online) with 1-200 traces highlighted over a base map. I'm happy to use online or offline tools (preferably Windows, Mac ok, last resort Linux) to achieve this, but the fewer steps and unfamiliar technologies the better. (I'm not really familiar with any offline GIS tools.)

Lots of other sites allow you to upload traces, which they'll display in a list - but I haven't found any that will actually render them all simultaneously.


texts: OSMNX KeyError when trying to compute edge_centrality on a graph composition

texts:
Context:

I am trying to merge several different graphs (some sharing nodes with others) downloaded from openstreetmap using omsmx.

While storing these different graphs in a list, I use the nx.compose_all() method from networkX to combine them into a single graph:

G_unique = nx.compose_all(Graphs_list)

Then I want to compute edge centrality on the graph edges using networkX's nx.closeness_centrality() function on its dual graph (by prior using nx.line_graph() function of the primal graph):

e_c = nx.closeness_centrality(nx.line_graph(G_unique))

All this is also well described by Geoff Boeing in its tutorial: 08-example-line-graph.ipynb on github.

The error:

But while trying to retrieve the value of that centrality, I always encounter the same error:

[e_c[edge + (0,)] for edge in G_unique.edges()]

Which returns a KeyError:

KeyError                                  Traceback (most recent call last)
<ipython-input-92-6410c27e0882> in <module>
----> 1 [e_c[edge + (0,)] for edge in G_unique.edges()]

<ipython-input-92-6410c27e0882> in <listcomp>(.0)
----> 1 [e_c[edge + (0,)] for edge in G_unique.edges()]

KeyError: (292494260, 244429342, 0)

As if this edge had never been found...

Precision: the same problem appears if I use nx.compose() in a for loop on all the graphs composing the list Graphs_list.

Question:

How can I solve this problem?

Environment:

Ubuntu 5.0.0-37-generic x86_64 GNU/Linux
Python 3.6.9 (default, Nov  7 2019, 10:44:02)
numpy 1.17.4
networkx 2.4
osmnx 0.10


texts: Optimizing SLD map styling implementation

texts: Styling a map on demand can be a great bottleneck if you don't have a good strategy to deal with it. I'm sure there is a lot of techniques to implement the styling of maps, but I would like to ask you if you know anything about the current state of research on this topic, and which are the most optimized "server" implementations of this.


texts:
I cannot open jp2-files with gdal:

This works fine (the Sentinel2 metadata):

~/copernicus/S2A_MSIL2A_20170111T100351_N0204_R122_T33UWP_20170111T100351.SAFE$ gdalinfo MTD_MSIL2A.xml

Driver: SENTINEL2/Sentinel 2
Files: MTD_MSIL2A.xml
Size is 512, 512
Coordinate System is `'
Metadata:
  AOT_RETRIEVAL_ACCURACY=0.0
  BARE_SOILS_PERCENTAGE=0.000388
  CLOUD_COVERAGE_ASSESSMENT=0.1018
  CLOUD_COVERAGE_PERCENTAGE=0.001909
  CLOUD_SHADOW_PERCENTAGE=0.389902
  DARK_FEATURES_PERCENTAGE=6.691170
[...]

But this not (one of the jp2 files):

~/copernicus/S2A_MSIL2A_20170111T100351_N0204_R122_T33UWP_20170111T100351.SAFE/GRANULE/L2A_T33UWP_A008127_20170111T100351/IMG_DATA/R10m$ gdalinfo L2A_T33UWP_20170111T100351_B02_10m.jp2
ERROR 4: `L2A_T33UWP_20170111T100351_B02_10m.jp2' not recognized as a supported file format.

gdalinfo failed - unable to open 'L2A_T33UWP_20170111T100351_B02_10m.jp2'.

I am using Anaconda 4.3.11 and created an environment with Python 2.7.13.
I want to open jp2-Files which are part of the Sentinel L2A-Product (converted from L1C with sen2cor).

I've installed openjpeg 2.1.0 and gdal 2.1.0 via "conda install":

> conda list

# packages in environment at /home/roland/anaconda2/envs/copernicus:
#
curl                      7.49.0                        1
gdal                      2.1.0                    py27_0
geos                      3.5.0                         0
geotiff                   1.4.1                         0
hdf4                      4.2.12                        0
hdf5                      1.8.17                        1
jbig                      2.1                           0
jpeg                      8d                            2
kealib                    1.4.6                         0
libgdal                   2.1.0                         0
libnetcdf                 4.4.1                         0
libpng                    1.6.27                        0
libtiff                   4.0.6                         2
mkl                       2017.0.1                      0
numpy                     1.12.0                   py27_0
openjpeg                  2.1.0                         5    sunpy
openssl                   1.0.2k                        0
pip                       9.0.1                    py27_1
proj4                     4.9.2                         0
python                    2.7.13                        0
readline                  6.2                           2
setuptools                27.2.0                   py27_0
sqlite                    3.13.0                        0
tk                        8.5.18                        0
wheel                     0.29.0                   py27_0
xerces-c                  3.1.4                         0
xz                        5.2.2                         1
zlib                      1.2.8                         3

The main problem seems to be that I don't manage to install the JPEG2000-driver:

> gdal-config --formats
gxf gtiff hfa aigrid aaigrid ceos ceos2 iso8211 xpm sdts raw dted mem jdem envisat elas fit vrt usgsdem l1b nitf bmp airsar rs2 ilwis rmf leveller sgi srtmhgt idrisi gsg ingr ers jaxapalsar dimap gff cosar pds adrg coasp tsx terragen blx msgn til r northwood saga xyz hf2 kmlsuperoverlay ctg e00grid zmap ngsgeoid iris map cals safe sentinel2 mrf wcs wms plmosaic wmts grib bsb netcdf kea hdf5 hdf4 gif jpeg png pcraster pcidsk rik ozi pdf arg

How can I do this?


texts: How to add heading to position using openlayers 3?

texts:
I've seen the geolocation with tracking example in OL3, what it does is it rotates the view depending on the heading of the device.

I've done some research on the web and all what I've found does the same thing

My question: is there a way to display the heading of the device without rotating the view ?


texts:
I have two polygon layers I'd like to overlay using QGIS or another open source tool.  What I want as output is a polygon layer with the extent of one polygon, but with each polygon in that layer subdivided into polygons based on overlapping polygons in the second layer.

This seems like a basic operation but I'm having trouble finding documentation on how to do it.


texts:
I use QGIS 1.4.0 to import a layer which is basically a road network converted from an .osm file to a PostGIS database table by osm2pgrouting.

I have skimmed through the documentation, but I am still not sure:

What happens when I edit features of such a layer in QGIS? Can I somehow write the changes back to the database?

Edit: This could be the answer, isn't it?


texts: How to add data slider bar control to Folium maps

texts:
I am making a choropleth map in Folium which displays data for certain age brackets. I would like to add a double ended slider letting users select the bottom and top of the age range to view. The data is from the census which has 14 different age brackets. Thus the slider should have discretized stopping points corresponding to the age ranges specified by the Census data. Something like this: https://digital-geography.com/filter-leaflet-maps-slider/
The only way I know how to do this in Folium would be to make each possible pair of ages its own layer and then switch between them in the layer menu. However this would increase the amount of data I would have by a factor of 98, and the map is already 270 KB with a single layer. So it would be an almost 20 MB webpage if I do it this way! Plus I would not have a slider.
Is there some other way of doing my desired functionality with an actual slider instead of the layer controls and without having to store such a large HTML page?


texts:
I need to apply the following query to 60+ tables:

UPDATE variable_table_name ab
SET    elem_nr = ef.elem_nr
FROM   ga_mit_elemnr ef
WHERE  (ab.elem_nr NOT IN (
       SELECT elem_nr FROM strassenelemente cd) OR
       ab.elem_nr IS NULL) AND ST_WITHIN(ab.geom, ef.geom)

Is there a way to automate this in PostgreSQL?


texts:
I'm new to Python, and am trying to install Qt/PyQt for PyQGIS. I downloaded the PyQt and SIP packages, and tried to install them but have run into errors.
Installing SIP: As far as I can tell, installation involves 3 command lines: 1. python configure.py 2. make 3. make install I managed the first two steps, but cannot figure why the third failed. I get the following error message:

Jocelynes-MacBook-Pro:sip-4.17 Jocelyne$ sudo make install
Password:
cp -f sip /System/Library/Frameworks/Python.framework/Versions/2.7/bin/sip
cp: /System/Library/Frameworks/Python.framework/Versions/2.7/bin/sip: Operation not permitted
make[1]: ** * [install] Error 1
make: *** [install] Error 2

I have no idea what Error 1 or Error 2 are, nor how to fix this problem.


texts: how to generate metadata file for semi automatic classification plug in?

texts:
I am working on satellite images. The region I am interested in is only a fraction of the whole satellite images. My laptop slows down tremendously when working on the whole image.
So I clipped the images.
Now I want to use the semi automatic classification plug in for classification but revive the following error message:

Error [8]: No metadata found inside the input directory (a .txt file whose name contains MTL)

I tried copying the original file which came with the original sized images when downloading into the folder of the clipped ones. This did not help.
How can i generate a new metadata file or work around this problem?


texts: Heatmap not Representing Distributions of Attributes Correctly

texts:
I have data for residential house sales in 2010 within England and Wales.

I've added x,y coordinates to the dataset enabling me to plot ~500,000 points across the country. From that I styled it into a heatmap which looks fantastic. Unfortunately the render always looks the same, despite changing the attribute by which points are weighed. Longitude is the only attribute where the map changes the weighting of points. See first 3 screenshots below

The task is to create a heatmap with points weighed by house sale price, not by abundance/concentration.

Interestingly if I use the interpolation plugin, lat (y) and long (x) attributes are the only ones which I can select. Same goes for the Heatmap plugin.

Finally, here is a list of my attributes. Is there an issue with the way I convert my excel work book into a tab deliminated text layer?


texts: Making list of features constructing dissolved feature?

texts:
I need to simplify my road network using the dissolve function to remove nodes not occurring on the intersections. the following code is doing that.

arcpy.Dissolve_management(self.address_cor_layer2, self.address_dissolved_ref_layer2, "", [["NAME", "FIRST"]], "SINGLE_PART", "DISSOLVE_LINES")

What I need is to keep track of features creating the dissolved feature. For that purpose I want to make a list of FIDs as one field in the dissolved feature.

Any idea how to do that?


texts: ArcGIS Data Interopabality Extension and Exporting from Mapcom M4

texts:
I am exporting from Mapcom's M4 database into Geodatabase or Shapefile.

I have the ArcGIS Data Interopability extension and I was able to connect and export the Lines, Points, polys, etc.

When I was reviewing the exports, there are quite a few features that still appear in Mapcom M4 that are not being exported. The Mapcom database is hosted on SQL  Server Spatial.

I really am at a loss why this might be happening. While identifying the features in Mapcom, I cannot see any discernible difference between the ones that are exporting and the ones that are not.

I am hoping that someone else has encountered this problem exporting from Mapcom, and if they have any ideas?


texts: R - Maps - How to add State border to City Map

texts:
I made a map of Brazilian cities using plot_brmap() and I want to show the state border.

But I don't know how to do it.

What I have been trying:

library(brazilmaps)
library(maptools)
library(RColorBrewer)

    # Define colors ----
    createColors  <- colorRampPalette(c("yellow","red","darkgreen"))

    # Define city and states (shapefiles) ----
    cidades_sp_mg <- get_brmap(geo = "City", geo.filter = list(State = c(31,35))
                               ,class = "SpatialPolygonsDataFrame")

    estados_sp_mg <- get_brmap(geo = "State", geo.filter = list(State = c(31,35))
                               ,class = "SpatialPolygonsDataFrame")

    # Create the city map ----
    load(file="dados/data.Rdata")

    dScore_data <- data %>% group_by(cod_municipio_ibge) %>% count()

    # Merge city and shapefile
    mapa_data <- plot_brmap(cidades_sp_mg,
                                       data_to_join = dScore_data,
                                       join_by = c("City" = "cod_municipio_ibge"),
                                       var = "n"
                                       )

    # Create state map
    mapa2_data <- plot_brmap(estados_sp_mg,
                                       data_to_join = dScore_data,
                                       join_by = c("State" = "code_state"),

    )

    MinhasCores_data  <- createColors(4)

    png("dados/Cities.png", width = 800, height = 600)

    print(mapa_data) +
      scale_fill_gradientn(colours=MinhasCores_data,name="Qtd.", na.value = "lightgray") +
      labs(title="Cities",caption = "Font") +
      theme(legend.title = element_text(face="bold"),title = element_text(face="bold"),
            plot.title = element_text(hjust = 0.5,size=20),
            plot.caption = element_text(hjust = 0.5,size = 10))

    dev.off()

Any suggestions ?


texts: Converting polygons to lines and find the intersection between the lines and polygons from another layer

texts:
I have three two different layers of polygons. I am using QGIS 3.16.2 (Hannover) in windows 10. I am not expert by using the vector data. However, I would like to find the intersection point between the polygons (Pink ones) and the buildings (Violet ones)
I was thinking first I can make the centerline of the seperated polygons (pink ones). Then I can extend the line to the polygon and find the intersection.
But I am kind of stoped in going through using QGIS. I would appreciate if you could give me some suggestion. I have tried to use the r.thin from grass, but I in QGIS I cannot create any output. I tried to define the center points of each seperated polygons (Pink ones) and connect them using the point to path. But since there is no order in the points which are created, it was also not possible.
I would happy if I get any suggestion either in QGIS or PyQGIS.
The question is if my overall thinking is fine (Make the centerline and extend to the polygon and then find the intersection)


texts:
I'm trying to write a Python script that imports arcgisscripting (stuck on 9.3.1 here) which is then run as CGI by a web app. It runs fine from the Windows command line but as soon as I drop it into the app it throws an Error 500: Internal Server Error. I checked the Apache logs and sure enough there's a ImportError: No module named arcgisscripting. All of my other modules are working fine (sys, os, etc.) Does anybody know what could be the issue here?

Update: I think it might have something to do with this: Why can't python find some modules when I'm running CGI scripts from the web?


texts: Dynamic oriented grid using QGIS 3 layout view (grid borders parallel to paper borders)

texts:
I use QGIS 3 to create maps with coordinate grid. Often, when I change the extent of a map, my coordinate grid tilts and becomes not parallel with the paper borders. This looks very awful, and I have to manually change the orientation of the grid.

Is there a way to automate this manual process and make the edges of my coordinate grid stay parallel to the edges of the paper?


texts:
I use a shapefile with CRS
("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0").
I want to calculate a buffer around the points with 'extract' from raster package. How can I know the units in this layer? Is it in meters? My idea is to calculate 1X1 km , thus buffer=1000. Is that correct?

try1 <- extract(raster,
                points,
                fun=mean, na.rm=T,
                sp = TRUE, # create spatial object
                buffer = 1000,
                method='simple',
                weights=T,
                normalizeWeights=F)


texts:
I have an issue withe the field calculator. It is impossible for me to multiplicate 2 columns (area * KZ). The result is always "0". I used the formula area * KZ....is this the right formula? LINK to shape file
Here the steps I did:

Using the "advanced python field calculator" to replace the strings "A"into 1, "B" into 2, "C" into 3 and so on
Using the field calculator (processing plugin) with the formula KZ * area without success


texts: ImageMosaic plugin for DEM rasters on Geoserver

texts:
I'm using 30m DEM rasters from the NRCS Data Gateway. I use image mosaic to merge them. But when i preview i see black gaps between rasters.

Coverager Parameters:

How can i remove that black gaps?


texts: Inserting text into field of attribute table in ArcGIS for Desktop?

texts: I am populating an attribute table. After adding a field i have to populate the field with different texts in a sequential order like water, air, car, school etc. i have completed inserting water using the field calculator, how do i add the others.


texts: Delete temporary files from geoprocessing in ArcGIS 10 and above

texts:
I am seeing my computers C: drive fill up quite quickly. I assume its from many geoprocessing things I have done. Is there a way to find the files Arcmap creates as temporary?

The folder I am in is, C:\Users\knappk\AppData\Local\Temp\


texts: How do you create a pie chart of a feature class showing the total area of each feature in it. I'm currently using ArcGIS Pro 2.6.1. When I click on the 'Create Chart' button in the Data tab for the feature layer, there is no option for a pie chart. In ArcMap, I would simply go to View>Graph>Create Graph and click on pie chart. Does ArcGIS Pro not have the function to create a pie chart outside of Symbology? I don't want every feature to show a pie, I just want one pie for the feature that I can move around the map as a graphic element.


texts:
I have ~100 million 2d vector topographic features (lines and polygons) in an Oracle 12c database.
I have 20 billion xyz points in an Oracle database table as number types - the points cover the same geographic area as the 2d vector data.
I want to height the 2d vector data based on the xyz point data - i.e. in practice this would probably involve creating a TIN based on the xyz points and then draping the 2d vector data over this.

Bearing in the mind the quantity of data involved, what technology options are there to perform this task?


texts:
I am new in GIS and I am trying to print a new map. If I use Google Street and web layer, it got printed without problem in the composer, but OSM Landscape does not appear on the printed document (it does on the screen though).

Any clue?

PS: I am using QGIS 2.4 on Ubuntu 14.04 64bits:
QGIS version    2.4.0-Chugiak   QGIS code revision  exported
Compiled against Qt 4.8.6   Running against Qt  4.8.6
Compiled against GDAL/OGR   1.10.1  Running against GDAL/OGR    1.10.1
Compiled against GEOS   3.4.2-CAPI-1.8.2    Running against GEOS    3.4.2-CAPI-1.8.2 r3921
PostgreSQL Client Version   9.3.4   SpatiaLite Version  4.1.1
QWT Version 5.2.3   PROJ.4 Version  480
QScintilla2 Version


texts: How to define the path which I surveyed 4 times?

texts:
I did a kinematic survey of a path. I automatically collected data with GPS every second while I was walking on the path (walking speed is not constant). I did 8 repetitions on the same path. Now I have 8 different point clouds for the same path. How can I define one point cloud from 8 point clouds?

The image shows how the points look in QGIS. I don't have enough knowledge to figure it out.


texts: Creating polygon regions with latitude and longitude values in qgis?

texts:
I've a set of places with their latitudes and longitudes, I've plotted them on QGIS but what I want is to create polygon from these points so that each point refers to a particular region around it. It's my first attempt at QGIS so I am completely new to it, if anyone could help?

Some of the points are :-

Dudhli  29.9713363  77.6049218
Raiwala 30.0254474  78.219132
Ajabpur 30.0431765  78.8255226
Fatehpur    30.0444012  77.7609435
Johdi   30.064216   77.5939097
Barontha    30.066753   79.0192997

I have around 106 such points.


texts: Problem transforming coordination system of this shape

texts:
I would like to convert this shape-file to WGS84 / UTM Zone 33N. But I failed in qgis.
What did I do wrong? Assigning the projection doesn't change it...
https://daswirdmanjawohlnochsagenduerfen.de/extern/airports.zip


texts: ArcGIS Experience Builder: Application does not show organizations basemap gallery

texts:
I am working with ArcGIS Experience Builder Developer Edition to build an application.
How I understand the docs and how it is described here, the basemap gallery shown in the application should be the one from the organisation (ArcGIS Enterprise).
However, this does not seem to work for me. Even though I have a basemap gallery defined in my Portal for ArcGIS (which works in the Portal), my application shows a different list (ESRI basemaps).
This is the list I get:

Any idea what I am doing wrong? Do I lack some configuration somewhere?
The connection between the application and my Portal seems to be fine, as I draw all my data and a Geocoding service from it successfully.


texts: Could not find function "coord_quickmap" in R

texts:
I was trying to actually crop part of a country (example is Korea) using this code in R:

korea + coord_quickmap(xlim = c(126.5, 128),  ylim = c(35, 37))

However, this error appears..

  Error in coord_quickmap(xlim = c(126.5, 128), ylim = c(35, 37)) :
  could not find function "coord_quickmap"

How can I do this?


texts: Extract manhole covers from ground LiDAR

texts:
I am wondering if anyone has come up with an automated or semi-automated workflow for extracting manhole covers from LiDAR data in ArcGIS? I was thinking about looking at the change in elevation relative to the road or perhaps the shape.

Any guidance?


texts: Passing input/output parameters to processing scripts using QGIS 3?

texts:
In the new qgis3 version it is still possible to add your own processing scripts to the processing toolbox. As I understand the discussion between the qgis3 developers (qgis3 is a very good work!!) the  way of passing input and output parameters has changed because the former way is not that pythonic, which is true. In qgis2.x input and output parameters were defined at the top of the script with a double # (comment) like:

##parameter1=file
##result=output raster
python code...

I've searched several hours in different forums, lists and so on, with no success how to start here and define input and output parameters...does anybody know where to get a documentation or example?


texts: Error installing/uninstalling GeoMedia WebMap Professional_6 on Windows 7 (64 bit)

texts:
On installing GeoMedia WebMap Professional_06.01.03.08 on Windows 7 Professional 64 bit. Following Error window pop up at the end of installation process or If I uninstall it.

I have .Net Framework 2, 3.5 and 4 installed.
Thanks!


texts: Deleting Duplicates from Multiple Tables in a GeoDatabase

texts:
I'm having trouble deleting duplicates from multiple tables in a geodatabase.  I have a field "facility_id" which should be unique, but I have multiple entries and need to delete all but one.  I've started a python script to delete the duplicates but I'm not sure how to finish it.  I tried a for loop to find all facility_ids in the values list and check for membership in the duplicate values list, but when I use cursor.deleteRow, nothing happens.

import arcpy, os
from arcpy import env

targetDB = r"C:\data\Temp\436e0250ef1748af92e7f1f75b5a0deb.gdb"
arcpy.env.workspace = targetDB
for table in arcpy.ListTables():
newTable = arcpy.Describe(table).basename
    if "ATTACH" not in newTable:
        with arcpy.da.UpdateCursor(newTable, "facility_id") as cursor:
            values = [r[0] for r in cursor]
            duplicates = [value for value in values if values.count(value) > 1]
            print "Original Table Values: " + str(newTable) + ": " + str(values)
            print "Duplicate Table Values: " + str(newTable) + ": " + str(duplicates)
            if values in duplicates:
                cursor.deleteRow()

output:

Original Table Values: ElecConduit: [u'PS145', u'PS022', u'PS022', u'PS108', u'PS134', u'PS999']

Duplicates Table Values: ElecConduit: [u'PS022', u'PS022']

Original Table Values: ElecLighting: [u'PS108', u'PS108', u'PS108']

Duplicates Table Values: ElecLighting: [u'PS108', u'PS108', u'PS108']


texts:
I added pyramidingPolicy into the batch function. However, it still uses
the median pyramidingPolicy. Something on the core batch process function must be needed/added?

batch.Download.ImageCollection.toAsset(rgb, 'to_TIMEMACHINE2_simple',
             {
             scale: 30,
             region:region,
             type: 'float',
             pyramidingPolicy: {'.default': 'mode'}
             })


texts:
I am attempting to run a GRASS command from the terminal and eventually a bash script, in which it links rasters in a folder to the Mapset previously defined. I was able to successfully create the Mapset running a single line of code, however I receive the following error when I link the rasters:

ERROR: Flag -e required also flag -c
Exiting...

It is to my understanding that I can run this code without having to create a new Mapset for each raster (reason for omitting -c flag), and can also run this without initializing GRASS 7.0 text (reason for including the -e flag).

My line of code that prompted the error is as follows:
grass70 /path/to/grassdata/urma_latest/PERMANENT/ -e --exec r.external input=urma_2018050900.tif output=urma1

My goal is to run this for a series of rasters housed in the folder /path/to/rasters/


texts: Get WMS features under a polygon in layer-vector

texts:
I have the buffer of a polygon in a layer vector, I want to know what features of a WMS layer are covered by the buffer, as the SLDSelect Control but with a polygon drawn.

Any ideas?


texts: Counting the occurrence of a specific letter within a string (QGIS)?

texts:
I have a field which contains a series of words separated by commas.

I am wondering if it is possible to count the number of times the letter 'G' occurs. Alternatively it might be possible to split the field by commas and go down another avenue...

Here's a few example fields:

Urban,Green,Urban,Green,Green
Urban,Urban,Urban,Urban,Green,Urban
Urban,Green,Green

Essentially I would like to produce a new column which tells me how many times 'Green' appears in each field.


texts: Unusual results for geoprocessing speed tests

texts:
I have been observing unusual performance with a Python geoprocessing script.  The (attached) script performs the following actions:

Use a search cursor to look up UTM zone corresponding to polygon
features
Create spatial reference object based on search cursor results
Convert .csv to feature layer and then to a point feature class

I've noticed markedly different processing times based on how the script is run:

32-bit processing using IDLE = 203 seconds
32-bit processing foreground script tool = 91 seconds
64-bit processing background script tool = 206 seconds

Why would this script perform so differently given the above conditions?  I certainly would not expect the 32-bit script tool run in the foreground to be 2X as fast as the other methods.

import arcpy, os, time

###IDLE Parameters
##fc = r'C:\path\to\polygon\fc\with\utm\zones\and\features'
##outws = r'C:\out\location'
##arcpy.env.workspace = r'C:\workspace'

####################
## Script tool parameters
fc = arcpy.GetParameterAsText(0)    # Feature class
outws = arcpy.GetParameterAsText(1) # Folder
arcpy.env.workspace = arcpy.GetParameterAsText(2)   # Workspace
####################

# Tables are .csv
tables = arcpy.ListTables()

start = time.clock()

# Look up which UTM zone .csv features are in
for t in tables:
    quad = t[7:17]
    print quad
    whereClause = """ "QUADID" LIKE '%s' """ % quad
    with arcpy.da.SearchCursor(fc, ("QUADID","ZONE"), whereClause) as cursor:
        for row in cursor:
            if row[0] == quad:
                utmZone = row[1]
                if utmZone == 10:
                    sr = arcpy.SpatialReference(26910)  # NAD_1983_UTM_Zone_10N
                elif utmZone == 11:
                    sr = arcpy.SpatialReference(26911)  # NAD_1983_UTM_Zone_11N
                elif utmZone == 12:
                    sr = arcpy.SpatialReference(26912)  # NAD_1983_UTM_Zone_12N
                elif utmZone == 13:
                    sr = arcpy.SpatialReference(26913)   # NAD_1983_UTM_Zone_13N
                else:
                    print "The UTM Zone is outside 10-13"
            else:
                pass

    # Convert .csv to feature class
    try:
        outLayer = "in_memory"
        # Now with the sr defined, create the XY Event Layer
        arcpy.MakeXYEventLayer_management(t, "x", "y", outLayer, sr, "z")
        arcpy.FeatureClassToFeatureClass_conversion(outLayer, outws, t[7:17])
        arcpy.Delete_management("in_memory")
        end = time.clock()
        print "In_memory method finished in %s seconds" % (end - start)

    except:
        # Print any error messages
        print arcpy.GetMessages(2)

print "Processing complete"


texts: Populate missing fields/labels for a grid of polygons with python?

texts:
I have an Esri polygon feature class in a grid that covers a geographic region, and one field, which we'll call 'Label' has an arbitrary label that consists of two digits for north/south, and two digits for east/west location. So when completely labeled a grid may look like this:

1010    1011    1012    1013
0910    0911    0912    0913
0810    0811    0811    0813

As you move east/west, the second pair of numbers increases/decreases respectively. As you move north/south, the first pair of numbers increases/decreases respectively.

The grids I'm dealing with may only have a handful populated, i.e.:

0000    1011    0000    1013
0000    0000    0000    0000
0000    0000    0811    0000   (with the zeroes just being null)

Or it may look like this:  (However, that it is never perfectly square like this.)

I've attempted to iterate through each polygon with arcpy, but that road I went down was too memory-intensive.

Can I save the 'Labels' field out into a Pandas or geopandas dataframe and manipulate it, and import back into the same field, keeping the geographic references/locations the same? I'm at a complete loss as to how to reference north/south/east/west (or a very literal up/down/left/right) sense of direction within a dataframe. How do I attack this problem


texts: Problem with Mapserver label encoding in Greek

texts:
I'm having a problem when trying to render a shapefile's labels that has Greek characters with the use of truetype and I'm getting weird characters. The code for my label section is as follows:

LABEL
                   TYPE truetype
                   FONT arial
                   SIZE 12
                   COLOR 50 50 50
                   PARTIALS FALSE
                   ENCODING ISO-8859-7
               END

I've also tried to use ENCODING GREEK  or even WINDOWS-1253 and UTF-8. But always getting weird letters.


texts:
I'm working on a bigger map for panaflex, 12'*10'. Since basemaps in ArcGIS cannot b exported as vectors (or may be I don't know how), I have to export a jpeg file or a pdf etc (pdf, even, without basemap has stripes of image). I have tried to export the file and I rec this error: "Invalid image parameter has been specified while rendering/writing the image". File gets exported in ai format but does not open in coral or adobe.

I would like to receive suggestions on how can I export this image to the best of quality with smaller size so that my designer may not have to work on this as exported map tends to be of huge size.

I do not receive any errors on page size A4 or A3 etc.


texts: Resample, ArcGIS: How do I avoid resample taking into acount NoData values?

texts:
I am trying to Resample a raster dataset to match the resolution of another dataset. I have fine resolution data on small islands (surrounded by NoData in the sea). When I re-sample the data to a resolution that is coarser than the island size, I lose this data. I need a way to re-sample to a coarse resolution, that takes into account the data on small islands, and not the NoData which takes up the majority of the new pixel size.

Any help?

David


texts: Geocoding to only reference Tileset using Mapbox GL JS?

texts:
I'm using Mapbox's Geocoding (search box) and was wondering if there's a way that it would only reference my Tileset data rather than the generic global address search?

Looks like Mapbox has a "Supplement geocoding search results from another data source", however the example is referencing geographic coordinates rather than addressing.

I'm still not sure if this is possible?

Basically I want to enter an address in the Geocoding search box and to have a particular polygon selected from my dataset on the map.


texts: How to connect to a SpatiaLite database using ogr in Python?

texts:
I need to connect to a SQLite/SpatiaLite database to copy data there. I have done the same with PostgreSQL/PostGIS like this:

connection_string = "dbname={} user={} password={} host={} port={}".format(
                self.dbname,
                config.get_config_value(dbsettings, "user"),
                config.get_config_value(dbsettings, "password"),
                config.get_config_value(dbsettings, "host"),
                config.get_config_value(dbsettings, "port")
            )
ogr.Open("PG:" + connection_string)

I followed instructions here where it is specified how to connect to PostgreSQL. But in the SQLite section, there isn't any information on connecting to the database. I suppose there also will be a connection string with a path to the database file instead of port and host but I need the exact syntax.

I am using Python 3.5.

EDIT:

I need to use ogr.open (not sqlite3.connect) because I store data to the database using CopyLayer (a method of ogr.open) as can be seen here (this is for PsotgreSQL):

    dsc_in = ogr.Open(file_name)
    dsc_out = ogr.Open("PG:" + self.target)

    layer = dsc_out.CopyLayer(dsc_in.GetLayer(), identifier,
                              ['OVERWRITE=YES'])
    if layer is None:
        raise Exception("Writing output data to the database failed.")


texts: Can pseudo-rectangular equal-area map projection have equidistant parallels?

texts:
Among the six pseudo-rectangular Eckert projections, three are equal area and three have equally spaced latitudes.  The equal area maps have parallels spaced to preserve area, I wonder if the same can be achieved by changing the shape of the meridians.

Is it possible for an equal area pseudo-rectangular map projection to be both equal area and have equally spaced latitudes?

If yes, what are some examples of projections that share all three properties (equal area, pseudo-rectangular, equally spaced parallels)?


texts: Merge features before creating new shapefile in PyQGIS

texts:
I'm working with QGIS 3.10 A Coruña on MacOS 10.13 environment.
I've got the following code (provided by @Kadir Şahbaz in this thread) to create new shapefiles from a loop selection of features for each different attribute value
layer = qgis.utils.iface.mapCanvas().currentLayer()

field_name = "CODNUT2"
idx = layer.fields().indexOf(field_name)
unique_values = layer.uniqueValues(idx)

folder = "/path/to/folder"

for value in unique_values:
    file_path = f"{folder}/_{value}.shp"
    layer.selectByExpression(f"{field_name}='{value}'")

    QgsVectorFileWriter.writeAsVectorFormat(layer,
                                            file_path,
                                            'utf-8',
                                            layer.crs(),
                                            "ESRI Shapefile",
                                            onlySelected=True
                                            )

Now I need to merge all features in each selection before creating the shapefile. I've tried including the following chunk of code in line 12
layer.combine(layer.selectByExpression(f"{field_name}='{value}'"))

obtaining the following error message
AttributeError: 'QgsVectorLayer' object has no attribute 'combine'

How can it be mended?


texts: How to create a raster of equidistant points with distance in meters using QGIS?

texts:
I need a set of equidistant points (like a raster) each in a 1km distance to the next over a certain area, in my case Germany. I know the feature "regular points" in QGIS where I can do just that but the distance you can enter there is not in meters but in coordinates or something.

I tried it in WGS 84 as well as in EPSG:31467 (Gauss-Kruger, common in Germany) which is supposed to be in meters but the results seem to be the same. I just can't find a way to set a certain distance in meters.


texts:
I am trying to create 3d buildings in ArcMap. I added SRTM, than created shape and draw on top of it few polygons and than specified to every polygon height (add and fill column height).

But when I am adding result to ArcScene polygons are drawing under DEM. But it should draw on it and 3D model should be on top.


texts:
I have a raster file with LCC projection  and a point dataset with SRID 4326. I want to extract raster value to point. When I use the query

SELECT ST_Value(rast, geom) val
FROM testraster5, pft_account
WHERE ST_Intersects(rast, geom)

I'm not getting any result. However the point lies on the raster. How do I solve this?


texts: Setting variable to search cursors to calculate percent using ArcPy?

texts:
I made a Python script tool and to finish my script I want to perform a percentage calculation using the 2 results from my search cursors. And the print the result to screen.

Here's what I have (but not working):

# Run search cursor on areas with slope code 3
with arcpy.da.SearchCursor(shapecalc_input, ['POLY_AREA'], "SlopeCode in (3)") as cur1:
     sum = 0
     for row in cur1:
       sum += row[0]

arcpy.AddMessage(sum)

# Run search cursor on areas with slope code 1 and 2
with arcpy.da.SearchCursor(shapecalc_input, ['POLY_AREA'], "SlopeCode in (1, 2)") as cur2:
     sum = 0
     for row in cur2:
       sum += row[0]

arcpy.AddMessage(sum)

arcpy.AddMessage((cur1/cur2)*100)

My search cursors give me the right values, just the calculation isn't working (last line).

Is my syntax wrong? or should I be using a different method completely? I'm new to Python.


texts: Drop down menu for Coded Domain not available in edit session

texts: Using ArcGIS Desktop 10.2.2, Advanced license.  I created a new coded domain for a string field. I added the domain to an existing field.  When I edit the features, i would expect to get a list of the coded values from the domain as a dropdown menu.  This is not happening.  What am I missing - is there a simple check mark or radio button setting?


texts: Create Memory layer in pyQGIS with projected coordinates

texts:
I am trying to make a memory layer in QGIS using python and then generate a grid base on calculations in meters. When I run the script, a banner comes and says the CRS is not defined. Then when I look at the derived information I see that the area of the Geometry and its position is wrong.

vector_layer = QgsVectorLayer('Polygon?crs:epsg=3879&field=ID:integer(3)&field=value:double','test',"memory")
poly = QgsGeometry.fromPolygonXY([[
    QgsPointXY(34.18639744, 95.15668787),
    QgsPointXY(12.0161471 , 82.35668787),
    QgsPointXY(24.8161471 , 60.18643753),
    QgsPointXY(46.98639744, 72.98643753)])
f = QgsFeature(vector_layer.fields())
f.setGeometry(poly)
vector_layer.dataProvider().addFeatures([f])
QgsProject.instance().addMapLayer(vector_layer)

This code will generate a square, but the location is wrong.
I tried to define CRS afterward in the code:

crs = vector_layer.crs()
crs.createFromId(3879)
vector_layer.setCrs(crs)
QgsMapLayerRegistry.instance().addMapLayer(vector_layer)

But it did not work either.
Can you please tell where I am making mistake?


texts: Geocoding sensitive addresses without API

texts: I have a large set of sensitive addresses that I need to Geocode on a server.  I cannot send them to an external API due to HIPAA restrictions.  Ideally I'd love to use an open source solution, even if it is not perfect.  Any ideas?


texts: Excluding multiple attributes from being labeled in ArcMap?

texts:
I am using ArcMap 10.1, and i am labeling a roadway system but would like to exclude multiple road labels from being displayed. I know how to exclude one road at a time using SQL Query by typing road<>4, but I can't figure out how to type that in order to exclude more than one road.

Any suggestions?


texts: Reclassify a raster value to -9999 and set it to the nodata value using python and or gdal

texts:
I modified a script which reclassifies a raster. It is based on the script by Roger in his blog (the second script) and it works pretty well.The problem I am facing is as follows:
 I would like to reclassify a value (in this case 100) to -9999 and then set -9999 as the nodata value. Below is the list of values to classify and the list to which the corresponding value should be reclassified. Note that value 100 should go to value -9999.

  classification_values = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,100] #The interval values to classify
      classification_output_values = [4,1,2,3,4,5,5,6,7,8,9,11,12,13,14,15,16,17,18,10,14,14,5,16,-9999] #The value assigned to each interval

However, when I view the resultant raster after running the reclass script, value 100 is not -9999, it defaults to 0, but -9999 is set as the nodata value.

In the script I added after line 33 in Roger's script

  dst_ds.GetRasterBand(1).SetNoDataValue(-9999)

If I reclass value 100 to 0 and set 0 as the nodata value the resultant raster is correct and everywhere were value 100 used to be there is nodata.

Not sure if the followong info might be of use:
Input raster is 8bit unsigned (have tried 16bit signed with no luck);
I have changed the numpy array from 8bit to 16bit signed (line 47 in Roger's script);
Output raster is 16bit signed.

How would I go about solving this problem?


texts: Help finding Postcode Boundaries for The UK, Spain, Canada, Ireland, Brazil?

texts:
I would be forever grateful for anyone's help finding current Postal Code shapefiles for any or all of these countries.

Spain
Brazil
UK
Ireland
Canada

Free is obviously better but because I may use them commercially I might need to pay for licensing.


texts: if-else & arcpy.Exists to check if a raster file exists and then execute code

texts:
In a geodatabase, I have several rasters with different suffixes such as (_a, _b, _c). Some gdbs contain all three rasters and others are missing one or more. I'm trying to use if-else statement along with arcpy.Exists to check if true then use it as an input, if false ignore it. Much of the script has map algebra equations and other raster tools. The issue is that with my current method, I'm unable to append the remaining strings from elif statements to the list so it can be used as the input for the mosaic to new raster tool. I'm using 10.4 and yes this is a crude method of doing of if-else statements and any suggestions are welcomed. This is my sample code block:

import arcpy,os

# Local variables:
rasA = os.path.join("C:/Test_Folder/","county.gdb/" + "raster_a")
rasB = os.path.join("C:/Test_Folder/","county.gdb/" + "raster_b")
rasC = os.path.join("C:/Test_Folder/","county.gdb/" + "raster_c")

print ("rasA", arcpy.Exists(rasA))
print ("rasB", arcpy.Exists(rasB))
print ("rasC", arcpy.Exists(rasC))

('rasA', True)
('rasB', True)
('rasC', True)

inrasList = [ ];
if arcpy.Exists(rasA):
    inrasList.append('rasA')
if arcpy.Exists(rasB):
    inrasList.append('rasB')
if arcpy.Exists(rasC):
    inrasList.append('rasC')

print (inrasList)
['Raster(rasA)','Raster(rasB)','Raster(rasC)']

# Step 1 - Mosaic to New Raster
rasLocation = "C:/Test_Folder/output/"
arcpy.MosaicToNewRaster_management(inrasList,rasLocation,"new_mos_ras","", "4_BIT", "1", "1", "LAST", "FIRST")

The append tool works using Brian's correction,  but the list is not recognized as a valid input for mosaic to new raster tool.

semiC = ";"
inrasList1 = semiC.join(inrasList)
rasLocation = "C:/output/"
arcpy.MosaicToNewRaster_management(inrasList1,rasLocation,"new_mos_ras","", "4_BIT", "1", "1", "LAST", "FIRST")

This successfully executes the tool, but creates a blank raster.

The code has been updated. I was using Raster(rasA) whereas the input should be simply rasA.


texts: I'm either missing something or failing at some basic math here, but is it possible to calculate the map resolution (map units/pixel) based on altitude (from sea level)? I have center coordinates and projection.


texts: Re-locating points randomly at certain distance in QGIS

texts:
I have GPS coordinates which identify aid projects in Africa. Now, I want to re-locate these projects, such that each project is re-located by 100 km in a random direction from its original location.
(Extra: I have country layers; can I force QGIS to re-locate only within the original country?)
I tried the Virtual Layer & Geometry Generator, both didn't work so far.
Any leads?
I use QGIS 3.14.15.


texts: Why do some tools / scripts not open for use in the ArcGIS?

texts:
I am using the Image Classification toolbar to classify a composite LANDSAT image and every time I go to the "classification" drop down and chose the "Iso Cluster Unsupervised Classification" option the tool fails to open and gives the following message:

"Failed to open tool IsoClusterUnsupervisedClassification (Spatial Analyst Tools)"

Also, when I try to run the tool using the search option and then choose the script for this tool that fails to open.

It is strange, every other tool within the tool bar works - I need this tool to give a raster output. I am using ArcGIS 10.1.


texts: How to properly reuse EE module more then once

texts:
From https://medium.com/google-earth/making-it-easier-to-reuse-code-with-earth-engine-script-modules-2e93f49abb13 I'm learning how to use modules.

So I can have mdl1:

var timeField = 'system:time_start';
exports.timeField = timeField;

Then on mdl2 I can:

var mdl1 = require('users/user/rs:mdl1');

Then I have mdl3 which required mdl2:

var mdl2 = require('users/user/rs:mdl2');

What if I need mdl1 on mdl3 as well, should I do the following on mdl3 again?

var mdl1 = require('users/user/rs:mdl1');

I put a print("module numnber") on each module and looks like the require call actually executes the script, I think by doing like above I'm executing mdl1 two times, which looks wrong. Isn't?


texts:
very early build ~ zip file

I am getting blank outputs for button images in a custom toolbar. Although the button information populates the dockable toolbar, they lack respective icons [see Fig.1]

Much of the following code was referenced from Add toolbar via Python

def initGui(self):
    # Add toolbar

    self.toolbar = self.iface.addToolBar("Graphic Groundwater")
    self.toolbar.setObjectName("Graphic Groundwater")

    # Create actions
    self.getDEM = QAction(QIcon(":/plugins/GraphicGroundwater/icons/getDEM.png"),
                           QCoreApplication.translate("Groundwater Modeling", "Get DEM Map"),
                           self.iface.mainWindow())
    self.getTopo = QAction(QIcon(":/plugins/GraphicGroundwater/icons/getTopo.png"),
                           QCoreApplication.translate("Groundwater Modeling", "Get Topo Map"),
                           self.iface.mainWindow())
    self.addBase = QAction(QIcon(":/plugins/GraphicGroundwater/icons/addBase.png"),
                           QCoreApplication.translate("Groundwater Modeling", "Add Base Area"),
                           self.iface.mainWindow())
    self.addNewBase = QAction(QIcon(":/plugins/GraphicGroundwater/icons/addNewBase.png"),
                           QCoreApplication.translate("Groundwater Modeling", "Add New Base Area"),
                           self.iface.mainWindow())
    self.addRiver = QAction(QIcon(":/plugins/GraphicGroundwater/icons/addRiver.png"),
                           QCoreApplication.translate("Groundwater Modeling", "Add River"),
                           self.iface.mainWindow())
    self.addStress = QAction(QIcon(":/plugins/GraphicGroundwater/icons/addStress.png"),
                           QCoreApplication.translate("Groundwater Modeling", "Add Stress"),
                           self.iface.mainWindow())
    self.addGrid = QAction(QIcon(":/plugins/GraphicGroundwater/icons/addGrid.png"),
                           QCoreApplication.translate("Groundwater Modeling", "Add Grid"),
                           self.iface.mainWindow())
    self.refreshGrid = QAction(QIcon(":/plugins/GraphicGroundwater/icons/refreshGrid.png"),
                           QCoreApplication.translate("Groundwater Modeling", "Refresh Attributes"),
                           self.iface.mainWindow())
    self.runModel = QAction(QIcon(":/plugins/GraphicGroundwater/icons/runModel.png"),
                           QCoreApplication.translate("Groundwater Modeling", "Run Model"),
                           self.iface.mainWindow())

    # Connect action signals to slots
    self.getDEM.triggered.connect(self.run)
    self.getTopo.triggered.connect(self.run)
    self.addBase.triggered.connect(self.run)
    self.addNewBase.triggered.connect(self.run)
    self.addRiver.triggered.connect(self.run)
    self.addStress.triggered.connect(self.run)
    self.addGrid.triggered.connect(self.run)
    self.refreshGrid.triggered.connect(self.run)
    self.runModel.triggered.connect(self.run)

    # Add actions to the toolbar
    self.toolbar.addAction(self.getDEM)
    self.toolbar.addAction(self.getTopo)
    self.toolbar.addAction(self.addBase)
    self.toolbar.addAction(self.addNewBase)
    self.toolbar.addAction(self.addRiver)
    self.toolbar.addAction(self.addStress)
    self.toolbar.addAction(self.addGrid)
    self.toolbar.addAction(self.refreshGrid)
    self.toolbar.addAction(self.runModel)

def unload(self):
    # remove toolbar on plugin unload
    del self.toolbar

[Fig.1]

The associated resource_rc file might also be of some diagnosis assistance...

<qresource prefix=":/plugins/graphicgroundwater/" >
    <file>GraphicGroundwater.png</file>
    <file>getDEM.png</file>
    <file>getTopo.png</file>
    <file>addBase.png</file>
    <file>addNewBase.png</file>
    <file>addRiver.png</file>
    <file>addStress.png</file>
    <file>addGrid.png</file>
    <file>refreshGrid.png</file>
    <file>runModel.png</file>
</qresource>


texts: GeoServer WPS request "Failed to retrieve value for input geom"

texts:
I am trying to do a buffer with WPS request builder of GeoServer.
I copied the URL of the gml3 output in GeoServer:

and I pasted it as the input URL

I ran the process and I received the following error:

  JTS:bufferBufferReturns a polygonal geometry representing the input geometry enlarged by a given distance around its exterior.Failed to retrieve value for input geom

the request is like this:

What can I do to fix it? Is the way I gave the URL as input of the process correct?

Should I change the name of the geom column of the input data?


texts: Automatically saving ModelBuilder output in the same folder it used for the input?

texts:
I am trying to create a model that takes an input layer, which will always be located in a unique folder on my hard drive, and spit out a set of analyses in the same unique folder the input was referenced from.

For example, I want to be able to add "X.shp" to the map from C:\Folder1, run the model, and have the outputs be saved into Folder1 automatically. If "X.shp" was taken from C:\Folder2, I would want the output to be saved automatically in Folder2.
Is it possible to automate this?


texts:
I want to do some advanced level scripting with VBA in the field calculator. I have done a simple coding in the code block. Here is the codes.

dim X
  if [Name]="University" then
    X=101101
  elseif [Name]= "College" then
    X=101102
  elseif [Name]= "Medical College" then
    X=101103
  elseif [Name]= "School" then
    X=101104
  else
    X=""
  end if

This is working correctly. But I want to do something more. Those are-

If the Name field contain the word University or College or Medical College or School then X return the expecting value.
Is it possible to create new fields through VBScript? If its possible then the X will split into two from its middle (3-3) and put the value in the newly generated column/fields.

If this possible it will help me a lot.

Thanks in advance


texts: How to calculate length of feature in MapInfo

texts: i Would like to calculate the Sum of pipe lengths in Mapinfo within a polygon. Any ideas?


texts: How to solve an installation problem requiring "liblwgeom" to work with "shp2pgsql"

texts:
I am trying to import shapefile to PostGIS. I have problem configuring shp2pgsql.

I am getting this error :

dyld: Library not loaded: @loader_path/../lib/liblwgeom-2.1.1.dylib
  Referenced from: /Library/PostgreSQL/9.3/bin/shp2pgsql
  Reason: image not found

I am trying to explore where the command is:

macbook$ locate shp2pgsql
/Library/PostgreSQL/9.3/PostGIS/man/shp2pgsql.1
/Library/PostgreSQL/9.3/bin/shp2pgsql

macbook$ ls /usr/bin/shp2pgsql
/usr/bin/shp2pgsql

I installed "liblwgeom" using homebrew,

macbook$ ls /usr/local/opt/liblwgeom/lib
liblwgeom-2.1.1.dylib

Do you have any idea how to resolve this issue .


texts:
I create a lot of county boundaries and was wondering if I could turn a column into a string so I can drop it into another statistical package. In this case SAS and R. For example, the CNTYPROPER field

would be returned as a string

("Sheriden","Mitchell")


texts: Unable to activate/toggle edit tool in QGIS 2.14

texts: I'm working with QGIS 2.14 Essen with Mac OS X 10.11.6 (El Capitan).  I was following along in the training manual where I was learning how to digitize a raster file (tiff format) provided by the training session.  I was able to enter the edit mode at one point, but now I find I can no longer activate (toggle) the edit tool.  I click on it, but no response. I've removed the raster file as a layer and tried to start over, but unfortunately, I had saved the project before and now I am stuck in this unactivated state. Any thoughts of how I get out of this state?  A screenshot is provided.


texts: Changing the label color with the Expression Dialog in QGIS 3.4.1

texts:
I want to separate two labels and then coloring each one separately, using QGIS 3.4.1.
So, I was trying to color a label within the Expression Dialog.
However, I could not figure out how to wright the function.

I tried all the Color function, but it did not worked.
The color's name was taken from the "Select Text Color" within the "layer Styling" bar and it a shade of green that supposed to replace the red.


texts: Snapping to vertex plane while editing in ArcMap?

texts:
Occasionally when I edit a new polygon feature I want to finish one side of it exactly where it meets the plane of an existing vertex of another feature (in order to create a 90 degree angle to that vertex).

This ability exists on AutoCAD (as can be seen in the picture attached), but I haven't found a way to do that on ArcGIS Desktop. I usually work around this problem by making bigger shapes than necessary and then trimming what I don't need.

Anybody know a solution?


texts: If/Then python script in Field Calculator not working

texts:
I'm working in Modelbuilder. I have a polygon in which I calculate the shape centroid in decimal coordinates (field called LONG_CENTR). Then I add a field (called FAJA) and I want to use a Field Calculator to assign it a value according to the centroid I've obtained previously. I think I have messed up the code in the field calculator.

Python-parser expression:

updateValue( !FAJA! )

Code Block:

def updateValue(FAJA) :
  if (70.500 <= match.fabs(LONG_CENTR) and match.fabs(LONG_CENTR) <= 73.500):
    return 1
  elif (67.500 <= math.fabs(LONG_CENTR) and match.fabs(LONG_CENTR) < 70.500):
    return 2
  elif (64.500 <= math.fabs(LONG_CENTR) and match.fabs(LONG_CENTR) < 67.500):       return 3
  elif (61.500 <= math.fabs(LONG_CENTR) and match.fabs(LONG_CENTR) < 64.500):
     return 4
  elif (58.500 <= math.fabs(LONG_CENTR) and match.fabs(LONG_CENTR) < 61.500):
    return 5
  elif (55.500 <= math.fabs(LONG_CENTR) and match.fabs(LONG_CENTR) < 58.500):
    return 6
  elif (53.500 <= math.fabs(LONG_CENTR) and match.fabs(LONG_CENTR) < 55.500):
    return 7
  else:
    return 0

The error I'm getting is:

ERROR 000539: Error running expression: replace( -69.256 )
Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 2, in replace
NameError: global name 'match' is not defined

Failed to execute (Assign Faja Value).

What is causing this error?


texts: Batch 'Create Feature Class > From XY Table' (ArcCatalog)

texts:
I'm using the 'Create Feature Class > From XY Table' tool in Arc Catalog, but have many data sets and would like to batch this process.

Is there an ArcMap Toolbox tool that will accomplish the same task? (Or a way to batch the task in Catalog?)


texts: Join a CSV file to GeoJSON using ogr2ogr in a Python script?

texts:
I have a CSV which I would like to join to a GeoJSON.

I found this answer which I tried to implement but I can't get it to work just yet.

Here's an extract of my python:

join_command = ['ogr2ogr', '-f', 'GeoJSON', '-sql',
                'SELECT wfp.*, wfp_desc.* FROM wfp.geojson left join "/home/wfp/wfp_desc.csv".wfpdesc ON wfp.vw_group_id = wfp_desc.vw_group_id',
                '/home/wfp/wfp_join.geojson', '/home/wfp/wfp.geojson']
subprocess.check_call(join_command)

But I get this error when I try to run it:

  ERROR 1: Unable to open secondary datasource `wfp' required by JOIN.

I've tried running it after removing the .geojson before 'left join':

join_command = ['ogr2ogr', '-f', 'GeoJSON', '-sql',
                'SELECT wfp.*, wfp_desc.* FROM wfp left join "/home/wfp/wfp_desc.csv".wfpdesc ON wfp.vw_group_id = wfp_desc.vw_group_id',
                '/home/wfp/wfp_join.geojson', '/home/wfp/wfp.geojson']
subprocess.check_call(join_command)

But then I get this error:

  ERROR 1: SELECT from table wfp failed, no such table/featureclass.

What am I missing?

Edit: 28/05/2019

Following the comment from klewis. I now have a working ogr command ready to go into python.

Using ogrinfo I found the correct layer name for the GeoJSON file. It was originally a shapefile and the original layer name carried over instead of the converted one.

However, following the join. The attribute fields have the layer names prefixed to them.

For example:
wfp.vw_group_id / wfp.name / wfp_desc.description

Is it possible to carry out the join without prefixing the layer name to the attribute fields?


texts:
I'm using QGIS 2.10.1.
Is there a way to avoid labels of different layers overlapping each other?


texts: QGIS getting OSM raster data to match with vector data

texts: I have downloaded OSM raster data from the openlayers plugin on QGIS, and put this on the map canvas using the web menu button.  The projection is EPSG 3395 WGS world Mercator.  However the projection for the vector data is EPSG 32633 WGS 84 UTM time zone 33N (Norway).  Using the properties button on each layer, I cannot find the options to get the same CRS for each layer.  The polygons in the vector layer seem out of place with the raster layer.


texts: Using shapely.affinity and affine_transform?

texts:
I have shapely.geometry Point:

from shapely.geometry import Point
point1 = Point(0, 100, 200)

and I want to swap coordinates Y and Z:

from shapely.affinity import affine_transform
matrix = [1,0,0, 0,0,1, 0,1,0, 0,0,0]
point2 = affine_transform(point1, matrix)

I expect point2 will be POINT Z (0 200 100), but I have POINT Z (0 0 100).

Is this a bug or my mistake?

python 3.7.0, shapely 1.6.4