texts: Making least-cost paths to avoid specific land-cover types using ArcGIS Desktop?

texts:
I am going to calculate least-cost paths toward a number of location, implementing Tobler's hiking function (i.e., considering slope as factor influencing speed of movement). I want to make sure that the paths will try to avoid some zones according to the land cover. The latter data is available as raster dataset.

Is it sound to reclassify the land cover types assigning a high value to the land-cover type I want the least-cost paths to avoid, and then multiplying the SLOPE by the reclassified LAND-COVER raster?

If it is ok, how should I proceed in case I want the paths to also avoid urban areas?


texts: What is the best tool to measure latitude the most accurately?

texts: I'm thinking it's the sextant but not sure. I tried looking it up on Google and came across this site but have found nothing useful. I am doing this for a research paper and need help. Thanks to anyone who may reply.


texts: Missing Polygons when importing a previously exported shapefile using QGIS (mmqgis plugin) Geometry Export to CSV.

texts:
I've obtained the following file from the US Census Bureau.

ftp://ftp2.census.gov/geo/tiger/TIGER2012/UNSD/tl_2012_42_unsd.zip

I've been working with these files and have noticed some issues with the results so I decided to investigate the tools I was using to see if they have something to do with the problem.

Steps to reproduce the problem:

Download the census shapefile from ftp://ftp2.census.gov/geo/tiger/TIGER2012/UNSD/tl_2012_42_unsd.zip.

Open shapefile with qgis.

Select Plugins>mmqgis>Transfer>Geometry Export to CSV file.

Pick a location to save the file.

Using QGIS, select Plugins>mmqgis>Transfer>Geometry Import FROM CSV file.

It will take a moment to run and then it will display the imported data.  The issue should be fairly evident as there are several large polygons that are missing.

Does anyone know what is causing this?  Thanks.


texts: Get location of PostGIS / GEOS topology exception

texts:
I use ST_LineMerge to combine many overlapping lines, splitting them where they branch off. A prerequisite to this process is to ST_Collect my geoms of interest, then ST_Node them to ensure all intersections are noded.

When running ST_Node(ST_Collect(geom)) on some of my geometries, I get the following error which I believe comes from GEOS here:

ERROR:  GEOSNode: TopologyException:
Iterated noding failed to converge after 6 iterations

When running this for large areas (all of the US, for example), it's difficult to identify where this is happening so I can fix the geometry. All input linestrings pass ST_IsValid(geom), and I end up having to narrow down potential locations by filtering my input and seeing what throws an error.

Is there any way to catch this exception and have it produce something like XXX problem at or near point....

Here's a MWE of my process:

with geoms as (
select  '0102000020110F000003000000DC874D65FCC25EC176032C6B350C5341B336429FFEC25EC1F962BBD9480C5341FC849518FFC25EC15BE20F5F500C5341'::geometry as geom
union
select '0102000020110F000006000000FA9BBFD3FCC25EC1B978232F390C5341B336429FFEC25EC1F962BBD9480C5341A77E6BE5FEC25EC1357C21334D0C5341C3EBA27BFEC25EC11BE5A4C34A0C5341B61D8CACFCC25EC1BCF273143C0C5341FA9BBFD3FCC25EC1B978232F390C5341'::geometry as geom
)

select  st_node(st_collect(geom)) from geoms;

Apologies for the ugly format- the process succeeds when using the ST_AsGeoJSON() versions of the polygons.

Here's my PostGIS full version:

POSTGIS="2.5.3 r17699" [EXTENSION]
PGSQL="110" GEOS="3.7.1-CAPI-1.11.1 27a5e771"
PROJ="Rel. 4.9.3, 15 August 2016"
GDAL="GDAL 2.2.3, released 2017/11/20"
LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.2.1" RASTER

How can I identify the location of this error (or feature IDs, etc) ?


texts: Import xls / xlsx / ods: field names are not in first row (just QGIS!)

texts:
I need to import a xls / xlsx / ods into QGIS. But while the first row (field names) of my XLSX in LibreOffice shows up on the first place as it should, QGIS shows the first row on place 452! As you can see in the screenshot below. I imported the XLSX to QGIS simply by drag'n'drop.

This issue seems to be very strange, because every program I tried, finds the row with field names on the first place! Also saving the XLSX to CSV results in the expected order: first row = field names (second screenshot). I tried saving and copying in all directions and variants - but QGIS still don't finds the field names on row 1.

Limitations:

I can not go for csv or any other format (just for testing)
the imported table need to be editable in QGIS.
The plugin "spreadsheet layers" is also not helping, because it uses virtual tables.
I can't share the data because there is a copyright.


texts:
PROBLEM: The functionality 'Insert your WMS/WMTS URL' in the 'Add a custom base-map' menu does not work and cannot be debugged properly.

BACKGROUND: Following the instructions at CARTODB MAKES IT EASY TO ADD THE PERFECT BASEMAP TO YOUR DATA VISUALIZATION I am attempting to add a WMS as a basemap. The WMS that I am using is dem_hillshades which I have copied from NOAA. Having tested the WMS in my browser (Chrome) and desktop GIS (QGIS) I can confirm that it works without problems in those clients.

ERROR: I am receiving the following error from carto.com 'The URL is either invalid or contains unsupported projections'.

First of all the URL is not invalid since it has been tested in the browser & desktop GIS and also having tried a GetMap request using the same link and adding the 'SRS=4326' parameter which is what carto.com supports does not address the problem.

Any ideas of similar experiences would be more than welcome.


texts: OpenLayers Multi-line drawing lines disappear

texts:
I've written code for plotting multiple linestring. I draw lines according to the distance entered and the number of lines. But in some cases, the lines start to disappear after a certain approximation to the line, but the main line never disappears. No operation or calculation and no errors. I couldn't figure out why.

The code example is as follows:

http://jsfiddle.net/OnurEk/ngv76ouL

What can be the reason for this and how can it be solved?


texts: element-wise subtraction of one property from another in a FeatureCollection

texts:
I would like to calculate RMSE for a model.

I have a FeatureCollection that contains the modelled and origianl properties. How can I do an element-wise subtraction (as the beginning of calculating RMSE)?

var model = ee.Classifier.cart().train(to_Train, 'N', bands);

var validated = to_Test.classify(model, 'N_');

var N = validated.select('N');
var N_ = validated.select('N_');
var sub = N.subtract(N_)

Produces

  N.subtract is not a function

This is what I would like to achieve

var RMSE = (N - N_).pow(2).mean().sqrt()

I think I have to make an array but can't seem to get there!

Here is a link to simplified version.


texts: Reducing number of connections required in FME

texts:
With FME 2017 I have to export multiple (205) feature types coming from an SDE Database. I have loaded them into Workbench which creates 205 separate Readers and I now have to create a writer for each one to export into a different format which is obviously very time consuming.

Which is the best method for outputting feature types individually as in this case?


texts: Converting/transforming coordinates from ETRS89 to WGS84 UTM 36N in QGIS

texts:
I am a student and just getting to know QGIS. At the moment I have to combine data which had been measured in ETRS89 with data in WGS84 UTM 36N. I created several vector layers for the data in ETRS89. Then I tried to transform the data into WGS84 by using "on the fly transformation" but it seemed to me that there was no change visible. The coordinates for one of our fixpoints are in WGS84 279000.000/4207000.000/1616.875 but QGIS shows me the coordinates 542000.00/4206000.00 (the height I can't even see).

How can I transform the coordinates?

I don't mean a change in projection I want to have a proper translation/conversion of the coordinates.

The more western points are the correct ones in the OpenStreetMap.  I think we still not have solved the problem. May be I should explain it a little bit better... I have got some data of Turkish students. They contain measurement points which had been acquired with a total station. We also measured with a total station in the same area but we measured in WGS84 UTM 36N relying on fixpoints which had been set witch GPS.

If you now look at the data, there has to be a difference because the coordinates of our Turkish colleagues begin with 542.../420... and ours 279../420... and also the height is divergent. But I'm not quite sure whether the information of the Turkish student was correct that it is ETRS89.

I really don't know how I can bring this data together.


texts:
I am pretty new to ArcGis Pro and would like to import (near surface geophysical data) which are available as x/y/z/p(arameter/value). The x/y locations are in UTM so no local coordinate system is needed to work around.
I am sure this is possible but couldn't find anything on the ArcGis ressources page.


texts: How can I use WFS layers with Configurable Map Viewer?

texts:
I have recently upgraded to v2.0 of CMV in the hope of using WFS layers. However, when configuring a WFS layer in viewer.js, no layers are loaded on the map (and there are no errors reported in the browser console).

Unfortunately I can't find any documentation around this functionality, and the only example I can find for loading WFS with the ArcGIS API for JavaScript doesn't seem to work.

Steps to reproduce:

Add a WFS layer to the operationalLayers array in viewer.js, for example:

{
            type: 'wfs',
            title: 'Roads',
            options: {
                url: "https://domain.com/arcgis/services/group/layer/MapServer/WFSServer" ,
                name: 'roads',
                version: "2.0.0",
                wkid: 29902
            },
            editorLayerInfos: {
                disableGeometryUpdate: false
            },
            legendLayerInfos: {
                exclude: false,
                layerInfo: {
                    title: 'Roads'
                }
            },
            layerControlLayerInfos: {
                exclude: false,
                noZoom: true
            }
        },

Load the map viewer.
Expected results:
I would expect to see the WFS layer on the map, as well as being able to switch it on/off it using the layer control widget.

Actual results:
The layer isn't added. Other operational layers do load, but the layer control is not visible in the left pane as it should be. Background layers also fail to load. There are no errors in the browser console.

I have also tried tweaking the code in _initLayer() within _MapMixin.js to load the WFS layer differently (as in the example here) but this does not load any features. It does execute a GetCapabilities request like the example, but I am not sure what to do next.

Environment:

Software    Version
CMV Version 2.0.0-beta2
Browser Chrome
Operating system    Windows 7


texts:
I have a task at hand to provide conflation functionality in our software.

We are using JTS for other WPS operations. Open source JCS library is available but we need it as a WPS service integrated in GeoServer. I found that GeoServer has support for JTS, but couldn't find if it has support for JCS by default in latest version (I am using GeoServer 2.2.3). In case it doesn't is it possible to add JCS into GeoServer?

I found that a new WPS process can be added to GeoServer, on similar lines I added the jcs jar (jcs-1.0.1.jar) to WEB-INF/lib in GeoServer and started my jboss instance which is hosting it. However, after restarting jboss, the WPS GetCapabilities is not listing jcs processes. I am sure i am missing something here but don't know where.

I have the JCS opensource code; will it have to be modified to make it compatible with hosting in GeoServer as in the link above or is there a way where GeoServer WPS can upload an api and host it?

I am new to GIS specially these GIS protocols part but trying to connect the loose ends.


texts:
QGIS has the ability to display GeoRSS by following these steps:

Click on "Add Vector Layer..." icon
Make sure the Data Source is "File" and enter the GeoRSS URL with "/viscurl/" prefix on the GeoRSS (e.g. /vsicurl/http://earthquake.usgs.gov/earthquakes/catalogs/eqs1day-M0.xml).

The problem I'm having is that the canvas is not updated when a new earthquake happens.  I have to quit qgis, restart qgis, then load the GeoRSS in order to see any updates.

Anyone knows how to show GeoRSS in real time?  I'd love to hear of any other FOSS solutions to this problem.


texts: Displaying points using Time property from geoJSON

texts:
Can I do something similiar to what https://openlayers.org/en/latest/examples/wms-time.html accomplishes here but with geoJSON as my layer source?

I imagine I should be able to still read in the whole geoJSON file but only display the points if they match my time variable and then continuously call map.renderSync().

I'm just lost on how to implement it, or if there's only an example of WMS because I need to be able to query my dataset.


texts: How to insert a statement (e.g. copyright) which remains visible in mapping window?

texts:
I have a QGIS workspace that I need to apply a copyright statement to, so we can show data to clients without breaching the terms of our data suppliers.

So far I've attached a label to a point in the bottom left of the area, and styled it nicely, but obviously this disappears when a user pans or zooms away from that point.

Is there a way to keep a block of text always visible in the QGIS window?


texts: Shapefile not displaying points (vertices) in ArcMap

texts:
this code finds all the vertices of every census tract for the state of California and Colorado. The code runs, and creates a shapefile with all proper filled fields (totpop, XCoord,YCoord, area1 etc.).
The problem now is that nothing is displaying when I add the shapefile on ArcMAP, however, when I open the attribute table, all the fields are filled.

fp = r"C:\Users\mr276170\Desktop\EXPLOSION!"

import os
import arcpy
import math
from arcpy.sa import *
from arcpy import env

env.workspace = fp

arcpy.env.overwriteOutput = True

#########
equal_shapeIN = fp + r'\CAandCO.shp'
out_shp = fp + r'\TheVertsAREAPOP.shp'
geometry_type = "POINT"
has_m = "DISABLED"
has_z = "DISABLED"

projection = arcpy.SpatialReference('USA Contiguous Albers Equal Area Conic USGS')

arcpy.Project_management(equal_shapeIN, out_shp, projection)

#########

arcpy.CopyFeatures_management(equal_shapeIN, out_shp)

arcpy.CreateFeatureclass_management(fp, 'TheVertsAREAPOP.shp',geometry_type, 'CAandCO.shp',has_m,has_z,projection)

dropFields = ["romney12", "obama12", "t_area", "totarea", "total_area"]

arcpy.DeleteField_management(out_shp,dropFields )

arcpy.DeleteRows_management(out_shp)

#########

fields = [('XCoord',"DOUBLE"), ("YCoord","DOUBLE"), ("area1","DOUBLE")]
for fields in fields:
    arcpy.AddField_management(*(out_shp,) +fields)

with arcpy.da.UpdateCursor(equal_shapeIN,["OID@", "SHAPE@", "STATE_NAME", "NAME", "totpop", "POLY_ID", "OBJECTID", "STATE_FIPS", "CNTY_FIPS", "FIPS", "FIPSnum","FIPS_NUMER","SHAPE@AREA"]) as cursor:
    #####add area

    totalPOP = 0
    for row in cursor:
        print("Feature {}:".format(row[0]))

        partnum = 0
        for part in row[1]:
            print("Part {}:".format(partnum))
            for pnt in part:
                if pnt:
                    Xcoord = ("{}".format(pnt.X))
                    Ycoord = ("{}".format(pnt.Y))
                    stateNAME = row[2]
                    countyname = row[3]
                    totalPOP = row[4]
                    idpoly = row[5]
                    idobject = row[6]
                    stateFIPS = row[7]
                    countyFIPS = row[8]
                    fips = row[9]
                    fipSnum = row[10]
                    fipsNumer = row[11]
                    totarea = row[12]

                    print Xcoord, Ycoord, stateNAME, countyname, totalPOP, idpoly, idobject, stateFIPS, countyFIPS, totarea

                    cursor1 = arcpy.da.InsertCursor(out_shp,["XCoord","YCoord","STATE_NAME", "NAME", "totpop", "POLY_ID", "OBJECTID", "STATE_FIPS", "CNTY_FIPS", "FIPS", "FIPSnum","FIPS_NUMER","area1"])
                               cursor1.insertRow([Xcoord,Ycoord,stateNAME,countyname,totalPOP, idpoly, idobject, stateFIPS, countyFIPS, fips, fipSnum, fipsNumer, totarea])

del cursor, cursor1

print "this works yo!"


texts:
I am using QGIS 3.12 and I want to update the content of a field by storing the feature's X and Y coordinates. The following code works given that a point layer is in the layer panel. The function works when the layer.updateFeature(feat.id()) statement is blocked, and so I could not understand the use of this method?
layer = iface.activeLayer()
layer.startEditing()
all_features = layer.getFeatures()
for feat in all_features:
    geom = feat.geometry()
    x = geom.asPoint().x()
    y = geom.asPoint().y()

    new_x = {feat.fieldNameIndex('x'): x}
    layer.dataProvider().changeAttributeValues({feat.id(): new_x })

    new_y = {feat.fieldNameIndex('y'): y}
    layer.dataProvider().changeAttributeValues({feat.id(): new_y })

    #layer.updateFeature(feat.id()) #?

# Commit changes
layer.commitChanges()


texts: QGIS issues on Mac OS Big Sur 11.2 GRASS and SCP issue

texts:
I have been having issues with QGIS on my Mac. The current operating system is Big Sur 11.2. The current workaround to the SCP plugin is to revert back to version 6 via these steps: https://fromgistors.blogspot.com/2020/10/SCP7.html
However, this method does not consistently work for others with older Mac OS like High Sierra. I am teaching a course and it doesn't always show up. Only version 7 shows up.
Per other recommendations:
I have tried uninstalling everything and reinstalling via Kyngchaos: https://www.kyngchaos.com/software/qgis/
This also includes the older release on the website like 3.10 as well as the newer releases. Each time, I uninstall everything including python and gdal and reinstall from scratch starting with python. The Processing plugin is not turned on and every time I try to activate it, QGIS crashes.
Downloading QGIS from the main website (https://www.qgis.org/en/site/forusers/download.html) and trying 3.16 and 3.18, I am able to get SCP 6 working with Big Sur but not with High Sierra.
I can't find older builds beyond 3.8 on QGIS website. They have them for Windows but not Mac OS for some reason.
I was trying to run operations instead in GRASS 7.8.4, but now that is no longer working. I am getting Unicode error and no amount of internet searching has helped me with this. I have no non-ascii characters in my path names that I can find so I am not sure what is going on.
Here is the error message regardless of version. To note, this is the i.pca run but this also happens with older tools like cluster:
QGIS version: 3.16.5-Hannover
QGIS code revision: 58ba7c1ed6
Qt version: 5.14.2
GDAL version: 3.2.1
GEOS version: 3.9.1-CAPI-1.14.2
PROJ version: Rel. 6.3.2, May 1st, 2020
Processing algorithm…
Algorithm 'i.pca' starting…
Input parameters:
{ '-f' : False, '-n' : False, 'GRASS_REGION_CELLSIZE_PARAMETER' : 0, 'GRASS_REGION_PARAMETER' : None, 'input' : ['/Users/Rachel/Documents/RS/Lab 7/ipca/band2after.tif','/Users/Rachel/Documents/RS/Lab 7/ipca/band3after.tif','/Users/Rachel/Documents/RS/Lab 7/ipca/band4after.tif','/Users/Rachel/Documents/RS/Lab 7/ipca/band5after.tif','/Users/Rachel/Documents/RS/Lab 7/ipca/band6after.tif','/Users/Rachel/Documents/RS/Lab 7/ipca/band7after.tif'], 'output' : '/Users/Rachel/Documents/RS/Lab 7/ipca', 'percent' : 99, 'rescale' : [0,2550] }

Traceback (most recent call last):
File "/Applications/QGIS-LTR.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Algorithm.py", line 423, in processAlgorithm
getattr(self, fullName)(parameters, context, feedback)
File "/Applications/QGIS-LTR.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Algorithm.py", line 499, in processInputs
self.loadRasterLayer(layerName, layer)
File "/Applications/QGIS-LTR.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Algorithm.py", line 746, in loadRasterLayer
self.setSessionProjectionFromLayer(layer)
File "/Applications/QGIS-LTR.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Algorithm.py", line 1044, in setSessionProjectionFromLayer
self.setSessionProjection(layer.crs())
File "/Applications/QGIS-LTR.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Algorithm.py", line 1051, in setSessionProjection
file_name = Grass7Utils.exportCrsWktToFile(crs)
File "/Applications/QGIS-LTR.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Utils.py", line 95, in exportCrsWktToFile
f.write(wkt)
UnicodeEncodeError: 'ascii' codec can't encode character '\xb0' in position 898: ordinal not in range(128)
Execution failed after 0.11 seconds

Loading resulting layers
Algorithm 'i.pca' finished
_________________________________

Does anyone have any ideas?


texts: Insert row into MapInfo table using c# and MapInfApplication object

texts:
I am attempting to insert a record, then eventually update a record using the MapInfoApplication object in c#.  I can currently read the table contents out into a listbox in a windows forms application no problem, however, writing is proving problematic.

MapInfoApplication myMI = new MapInfoApplication();
myMI.Do(@"open table ""\\servername\\directory path\\TestTable2.tab"" as testtable2");
myMI.Do("select * from testtable2 into selection");
int count = int.Parse(myMI.Eval("tableinfo(selection,8)"));
myMI.Do("fetch first from selection");
for (int i = 0; i < count; i++)
{
    comboBox1.Items.Add(myMI.Eval("selection.col1"));
    myMI.Do("fetch next from selection");
}

the above works fine, below does not

myMI.Do(@"insert into testtable2(picname) values(""pic 5"")");
myMI.Do(@"commit table testtable2");
//mostly get error like table being edited by another

Can anyone assist?


texts: QGIS 2.0: Automating tasks, text files and shapefiles

texts:
Wondering if it's possible to automate some repetitive tasks in QGIS (2.0).
This is my general workflow for importing some xyz type text files from another software package which export (not overly efficient) fixed column data:

import fixed column text file into Excel
add fieldnames
create unique id field and populate using a combination (concatenate) of two fields in the existing data
remove null rows in between data groupings
remove unwanted columns
save to csv
import into QGIS
convert to polylines and save as shapefile

I'm thinking that this is achievable in QGIS alone via some scripting?
Example File:

So, I want to keep the first 4 columns (X,Y,ID,Value), delete the final column, remove the rows with text, and create a polyline shapefile based on a field which is a combination of the 3rd (ID) and 4th (Value) columns (e.g. 0_2000, 0_2100, 1_2100, etc)
.


texts: Deleting more than one node at once in QGIS

texts:
I would like to delete some nodes in my layer. I've got them about 50 or even more.
I found some solutions here:

https://www.youtube.com/watch?v=lbgmNXXJCzA&ab_channel=JustinAdler
How can i delete multiple nodes, for a polygon, with Quantum GIS?

but they don't work.
Is there some option to make i.e. rectangular selection within the selected item and remove several nodes simultaneously when hitting the Del button on the keyboard?
I would like something like here: Remove a node in QGIS but for more nodes at once.


texts: org.geotools.jdbc.JDBCDataStore change column name from the_geom to shape?

texts: org.geotools.jdbc.JDBCDataStore stores the geometry for a shape file in a table column called the_geom. Is it possible to change the name of this column to something other than the_geom, perhaps shape for example?


texts: Updating set of attributes on vector layer from .csv table?

texts:
I have over 1,000 records in a .shp file, all with a unique ID field (LINKID), along with many other fields containing data (TXTFLD1, TXTFLD2, INTFLD1, INTFLD2, etc).

These records have been copied into a spreadsheet, checked, updated, and re-imported into QGIS as a .csv. The .csv file still has the unique LINKID to create a link between the two datasets, and they still have exactly the same field names.

I am fairly new to QGIS and I am completely new to Python.

How do I use it to update the existing vector layer attributes from the updated .csv file?

Or is there a better way of doing it?

I want to overwrite the data fields in the .shp file with the data in the .csv.
This is a process that I will have to do hundreds of not thousands of times, and I do not particularly want to have to update each field manually every time, and I would like the process to be automated.


texts: Digitizing Menu unavailable in QGIS 2.18.15

texts: I installed the latest version of QGIS 2.18.15 and the Digitizing menu seems to have disappeared in the properties options of vector layers. I want to enable the 'Use Last Attribute" feature


texts: Openlayers3 map stretch when resizing browser window

texts:
When resizing the browser window, the map image is enlarged or stretched in the resize direction when renderer is set to webgl or canvas, but not with DOM.
This can be seen on this example:
http://openlayers.org/en/v3.10.1/examples/side-by-side.html?q=render

Where does this comes from, and is there a way to prevent this appart with DOM renderer?


texts: Multiple calculations with if then statements and NULL values

texts:
I'm trying to calculate a total demand in ArcGIS 10.2.2, stored in fields Demand1, Demand2, ... , Demand20. Some of them contain a numeric value. The others have a <Null> value. The sum therefore always returns <Null>.

I've tried in Python to define, in the pre-script a function:

def D(a):
   0 if a else 1

(as recommended in Field calculator if and elseif statement with <Null> values in ArcGis 10.1, with correct indents, etc.)

then in the code

 D( !Demand1! ) + D( !Demand2! ) + ... + D( !Demand3! )

It cannot compute. Do you have any suggestions?


texts:
How can I change the offset of the CartoDB tooltip from the cursor? When I add a  to the template content, the cursor overruns the tooltip and I would like to tweak the offset a bit so its 5,-5 pixels up and left of the cursor.

vis.addOverlay({
   layer: spillLayer,
   type: 'tooltip',
   template: '
{{type}} spill of
{{totvol}} gallons
'
});


texts: Map Package or Relative Paths for ArcMap MXD export

texts:
Is there a way to use Map Packaging so that multiple mxds draw data from the same geodatabase when unpacked?  We've broken our project area into 14 divisions.  Each division has it's own geodatabase and there are about 30 mxds per division drawing data from that geodatabase.  I've been asked to give all our project data (over 400 mxds) to clients and I want to do that efficiently from a file storage standpoint?

In my limited experience with map packaging, it looks like it would create one geodatabase file for every MXD. If true, I'd have over 400 .MPK with each one having it's own geodatabase. Is there a way to reduce the redundancy in geodatabases since there are about 30 mxds pulling data from the same .gdb?

If file duplication and storage is a concern, would it be better to abandon map packages and provide the mxds and data using relative paths?


texts: Different area values from Zonal Geometry and Tabulate area

texts:
I am working on Landcover map for Ghana. I calculated the area of each landcover type with Zonal Geometry in ArcMap. Then I wanted to know the area of landcover in each region (as in location). But when I used tabulate area in zonal statistics, the sum of total area for each landcover type among the regions was different from the one I had in using zonal geometry.

Please what could be the cause?


texts: How to visualise "distance from amenities"?

texts:
I'm interested in calculating for every point the distance to the nearest amenity of each of a certain number of types (train station, supermarket and so on). I think the calculation will be straightforward, but I'm not sure how to visualise the result. A heatmap doesn't seem quite right.

Things I want the visualisation to do:

No distinction needed between a point that is near 1 supermarket vs near 2 supermarkets
Show each attribute individually (for maybe 3 attributes)
Use a step function for distances (eg, 0-499m is one colour, 500-999m is another...)

Is there any standard kind of visualisation for this kind of analysis? My stack is Leaflet, TileMill, PostGIS, OpenStreetMap, comfortable with Python tools but would consider branching out.


texts: Converting geotiff from geographic coordinates to UTM_WGS84 using Python?

texts:
Does anyone know how to quickly convert a number of geotiff images with geographic coordinate systems (GCS_WGS84) to UTM_WGS84? Preferably using python.

I tried gdalwarp below but it failed:

The spatial reference of the source image is shown below:


texts: Credential Dialog, disabling default Realm information messages disclosing password

texts:
I'm using the QgsCredentialDialog for determining user access to the database. When a user is failing to log in correctly the log in credentials are shown in the dialog (Fig. 1).

(success, uname, passwd) = \
 QgsCredentialDialog.instance().get(connInfo, None, None, "Log in")

Despite the credentials being faulty there's the case where username is misspelled but the password is correct.
As a minimum, I'd like it if the entered password was not shown in the "Realm" message at the top of the dialog.

Is it possible, via Python, to disable this messages in the dialog or am I better off creating my own credentials dialog in Qt?


texts: One widget hijacking another during on change event in Web AppBuilder app

texts:
I have 2 widgets: a Manage Folders widget and a Task Manager widget. They sit next to each other in the navbar of my Web AppBuilder version 1.2 app. There are several other widgets as well. If I launch Task Manager after any of the other widgets, there are no problems. However, if I launch the Manage Folders Widget, close it, and then open the Task Manager widget, I immediately notice that a dropdown is sized a little larger than normal. This dropdown controls what task items are displayed in a grid below it. Normally, if I select a new item from the dropdown, it triggers and on change event that loads the tasks for the new selection into the grid. Now however, when I select a new item from the dropdown, nothing happens; the items from the previous selection (or default) remain the same.In investigating this issue in Chrome DevTools, I noticed that when I make a change in the Task Managers dropdown, the on change event gets triggered but then it goes to the 'loadData' function in the widget.js file for the previously closed Manager Folders widget. Note: they both have a function called 'loadData'. How could this be?! What seems to be happening is that the "thisWidget" variable is still assigned to "ManageFolders" so that when I call:

dijit.byId('subprojectSelect').on('change', function (e) {
    thisWidget.loadData(e, proj_type_obj);
});

in TaskManager, it is going to the "loadData" function in ManageFolders. How can I kill the ManageFolder "thisWidget" and make sure that it is assigned to TaskManager instead?

Any suggestions (other than renaming one of the loadData functions? Why would this be happening? Could this be something set incorrectly in one of the .json files?


texts: Convert raster cells to vector (polygon) grid?

texts:
I would like to convert a raster (representing a study area) to a polygon layer with features representing individual raster cells.

The raster layer contains roughly 30 million cells @150x150m.

I have tried using the Create Fishnet tool, but I want the cells to match the study boundary, not the study area boundaries EXTENT (i.e. many more unnecessary grid cells), and selecting the grid cells from the fishnet that intersect with the boundary area does not work (freezes) because of the large number of grid cells (~30 million).

Could anyone recommend a procedure for converting the raster to vector grid cells?


texts:
I have installed Tilemill to export my map as mbtiles,and I used mbutil to export as PNG images & UTFGrid.So forth,I have a Tilemill project with MSS(Carto) files.Now the data source is going to change weekly and I wish to automate this action.

So does Tilemill support command-line scripting or something else to export my project with given properties ??

Current System :
Windows 7 x64

EDIT in response to @Sunil:

  TileMill doesn't have any features for automatically manipulating your Carto stylesheets, but you could use a bash script to manipulate the Carto files directly. TileMill also has an export command that can be called from bash.

On Mac OS X, you could run a command like this in the Terminal for more details about using the export command:

/Applications/TileMill.app/Contents/Resources/node /Applications/TileMill.app/Contents/Resources/index.js export --help

The exact command depends on how you installed TileMill

Unfortunately,same command for Mac OS X does not work on my current system.What does the last sentence mean,is there different installations of Tilemill ?


texts: listing multiple files for a merge (arc)python

texts:
I have created a script to convert GPX into features with python. This script works but now I want to add all these separate feature classes (point features) in one featureclass/layer and convert it into a line feature. The following is the coding I did to merge the layers.

import arcpy
import os

arcpy.env.overwriteOutput = True
MergeFolder = "D:\Altwym_2017\programmeren\Geo_data.gdb"
Total = "D:\Altwym_2017\programmeren\output\total"
arcpy.env.workspace =MergeFolder

Mylist= arcpy.ListFeatureClasses(feature_type= "point")

arcpy.Merge_management (Mylist, Total)

The following Error occurs;

  Traceback (most recent call last):   File
  "D:/Altwym_2017/programmeren/scripts/Merge.py", line 13, in
      arcpy.Merge_management (Mylist, Total)   File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\management.py", line 4124, in
  Merge
      raise e ExecuteError: ERROR 999999: Error executing function. Failed to execute (Merge).

I think the listing is the problem since no error occurs if I replace the Mylist with the actual filenames of the converted GPX files. Does anyone know how I should change my way of listing (or use some other function/command to select a huge amount of files).

Filling in/ using some function for finding the right names is not possible since the names have nothing in common.


texts:
I have a table with a column named "Kategorie", filled with values. I added a new column named "RULE". I want to fill this column with Strings depending on "KATEGORIE". I have a PY-Script and it doesn't work.

What is wrong?

Expression:

"getKat(!KATEGORIE!)"

Codeblock:

"""def getKat(KATEGORIE):
    if KATEGORIE = "0":
        myval = "Autobahn"
    if KATEGORIE = "1":
        myval = "Bundesstrasse"
    else:
        myval= 3"""
arcpy.CalculateField_management(ERG, "RULE", expression, "", codeblock)

Both answers so far don't work.
I think the problem is the line "expression" or "codeblock". I try to print the result of "myval".The result is a error message "name myval is not defined"

expression:

"getKat(!KATEGORIE!)"

codeblock:

"""def getKat(KATEGORIE):
    if KATEGORIE == "0":
        myval == "Autobahn"
    elif KATEGORIE == "1":
        myval == "Bundesstrasse"
    else:
        myval == "3"
    return myval"""
print "myval", myval
arcpy.CalculateField_management(ERG, "RULE", expression, "PYTHON", codeblock)


texts: Gisgraphy results appear in two languages

texts:
I hosted Gisgraphy server for geolocalization issues and imported the data for Saudi Arabia. I need it mainly for reverse geocoding issues. When I request some coordinate it returns the result in two languages (English and Arabic).

I need the result to be in the local language (i.e. Arabic) or at least the city name should be into Arabic. How to achieve this?

edit:
example of a request and an answer :
to request I use the link:

http://localhost:8080/reversegeocoding/search?format=xml&lat=24.698246&lng=46.846404%22;

and the answer :


texts:
I am finding difficulty when it comes to retrieving date values in hive. I am using esri git hub tools for hadoop and my final objective is to display these columns in arcmap.
My query is

> create external table test1(DISPLAYSCALE int, CREATED_DATE date, LAST_EDITED_DATE date)
> ROW FORMAT SERDE 'com.esri.hadoop.hive.serde.JsonSerde'
> STORED AS INPUTFORMAT 'com.esri.json.hadoop.UnenclosedJsonInputFormat'
> OUTPUTFORMAT 'org.apache.hadoop.hive.ql.io.HiveIgnoreKeyTextOutputFormat';

When I try to use the select * query, I get the following error

  select * from test1 limit5;
  Failed with exception java.io.IOException:java.lang.ClassCastException: org.apache.hadoop.hive.serde2.io.DateWritable cannot be cast to org.apache.hadoop.io.Text

As per the json the datatype for  CREATED_DATE and  CREATED_DATE are esriFieldTypeDate and the values are in this format say 2013-11-20 09:39:25.000001.
So i used the date datatype while creating the table, copied it to HDFS using the unenclosed json and used the select * query to retrieve the columns, but I get the above error.
To get the values we are creating the same table with string data type respectively instead of date and we are able to get the values .
I tried it with different columns including shape attributes but the initial step of selecting the query itself is not working and am not able to proceed further

Can you suggest a solution for this problem. This question may seem silly but I am pretty new to programming.


texts: Finding the pixel coordinates using GDAL/numpy

texts:
Is there a way to get the XY cordinates of a pixel based on the array index using numpy/gdal?

I have a new set of extents for a raster and, based on these newer extents, want to determine the centroid of the raster (polygon). Wondering if there are any functions which can return the coordinate values based on the pixel indices.


texts:
I need Earth radius in pixels for use in Mercator projection. Well let's say that our map's height is 1478 pixels and distance between 0 pixel and 1478 pixel is 2.52 degree. Radius is distance from the equator to pole so the formula for count Earth radius in pixels is R=90/2.52*1478 which is 52785.71429. Well then pixel number for 48.54 degree is 51287.99461 and for 51.06 is pixel number 54885.98194 counted by Mercator so the diference is 3597.987328
 which is not correct because it should be those 1478 pixels. I think that Mercator formulas are good because when i make inverses of pixel i get the degrees as in input so mistake can be only in Earth radius. Am I missed something, are my calculations right?  Thanks in advance.


texts:
I have three type of layer and for each layer a different marker, which I styled in CSS. The problem now is, that if someone uses VoiceOver or something like this, it's reading the 'html'. In my case, just the first letter. Is it somehow possible to add the ARIA-label attribute to my marker, that VoiceOver will read feature.properties.name and not the html?
     pointToLayer: function (feature, latlng) {
            return L.marker(latlng, {
                icon: L.divIcon({
                    className: feature.properties.type,
                    html: feature.properties.type[0].toUpperCase(),
                    iconSize: L.point(50, 50)
                })  })    }

I butchered a bit in CSS:
 divm {
        color:transparent;
        font-size:0px;

    }

And used it like this in the html:
html: feature.properties.type[0].toUpperCase() + '<divm>' + feature.properties.name + '</divm>',

It's not the ARIA label, but VoiceOver is reading the name. But is there a more elegant way with ARIA?


texts: How to balance/equalize color of geotiffs before creating a mosaic

texts: I have a bunch of geotiffs and I want to balance/equalize the color between them, so when I create the final mosaic there will not be any "too obvious" artefacts. GDAL helps to create the mosaic (as shown in several posts) but how can I do the trick for the color? Any help is much appreciated.


texts:
I am working on a thematic site map for a scientific study. In addition to my own data that will be displayed on the map, I have included reference data from U.S. federal agencies (e.g. National Hydrography Dataset (NHD) by United States Geological Survey (USGS)). Because NHD is in the public domain, my understanding is that I am free to reproduce it, but I still must acknowledge USGS.

So, do I need to cite a public domain source? If so, what information must the citation include?


texts: Large survey on unconnected underground networks

texts:
We (read I) are currently surveying all the stormwater pits and pipes in our shire, and have run into a few issue with handling unknown network connections.   The problem is that a lot of our network is not pit-to-pit like any good network should be (I can thank lack of standards for that) so we have a lot of pits that have a pipe running out to another pipe somewhere underground and it’s just busted in with no junction pit.

We are able to get a invert on the outgoing pipe from the pit but not where it joins the main pipe, because well, we don’t know where that main pipe is or even how deep it is.   I think a picture will help with what I am describing:

So the green dot is the new surveyed pit and the X is the invert position.  We have a rough direction for the outgoing pipe but no idea how far back into the footpath the main pipe is, so creating a connected network is very hard.

Now imagine over half your network is like this, and this is an easy case.  How do you handle it?

I’m interested to know how other people would/have handled doing a full survey of their stormwater (or any other underground) network.  Our ultimate goal is to have a complete connected network with inverts at all joins in order to be able to run models.

I'm very interested to know how you would get inverts for the pipes that have no junction that comes to the surface.  Is digging each junction up the only option?


texts:
I need to create a set of points with a fixed spacing and given distance from a  a georeferenced coastline. My idea is trying to obtain a parallel line and then placing points on it.
To do the first I tried using the Offset lines tool, but the processing gets stuck/really slow: apparently it's because of the small distance I need to set and the complex geometry of the coastline.

Is there any faster way to do this?


texts: Error: Grid layer was not correctly generated

texts:
Faced this error already several times when I was playing with processing algorithms:

I used hi-res DSM/DEMs out of UAV surveys.

Resolution between 2,5 - 5 cm GSD. Output format: .tif

Any ideas how to fix it?


texts: I have a large polygon dataset (26,000) records with numerous attributes columns.  I wish to calculate the percentile rank of each attribute.  Say I have population_density in an attribute, is there a simple formula to assign a rank in another attribute field on the percentile score of the value in the first attribute column (pop density)?


texts: I would like to color code floors (like the one in the picture). How do I go about doing it? My 3d model is in glTF format. I am using cesium as webclient platform.


texts: I have Geoserver serving cached layers (WMTS) as PNG and PNG8. I am able to render my layers using the runtime in PNG without an issue, but every time I try and render in PNG8 I just get a blank screen despite serving an existing cache for both PNG and PNG8. Does anyone have any ideas as to how I can work around this? I am looking an expanding my layers and want to keep the rendering as light as possible, thus why I want to use PNG8.


texts: Using OTB-Wrapping of Orfeo Toolbox (Java)?

texts:
I can not find the code/binding/whatever for Java...

http://orfeo-toolbox.org/otb/

https://en.wikipedia.org/wiki/Orfeo_toolbox

Can someone help me on that issue?


texts: Truncated features in a feature layer(hosted) still consumes database storage space in ArcGIS Online

texts: I have a feature layer hosted in AGOL, it had 333 records and used 18MB of DB storage in ArcGIS Online. I then removed all the features using Truncate. It still shows me the same storage space even there are no features.


texts: How to detect duplicates and keep only one?

texts:
I've got a table with lots of points (200.000). Some of them are duplicates (not exactly but e.g. a metro station has eight entrances) and I want to keep only one of these metro entrances.

I am able to identify the duplicates with the following statement (with desc_string = Metro Station Name):

--table of just duplicates
create table dups as
select t1.* from ciss_poi as t1 join ciss_poi as t2
on st_dwithin(t1.geom,t2.geom,(0.5/111.111))
where t1.desc_string=t2.desc_string and t1.gid != t2.gid

Then I this:

The picture shows two metro stations with total 16 points.
My question now is how can I only keep one point of each metro station and delete the others?
Maybe creating bounding boxes around each of the doubles would be a hint, but I do not know how to do it.

EDIT: Maybe this might another way of solution:
Maybe one possibility is to draw bounding boxes around points where this point have the same name and that are in a maximum distance of 500 meters to each other.
I guess I can do with st_envelope but I do not know how to write the SQL query.


texts: QGIS3: Import csv_file with DateTime Values


texts: Show GML on map (using an example GML file)

texts:
How can I show features (from an example file)? I tried many, many ways. This is a different route using a standard GML example file from OpenLayers.
Problem : With the shown code example, I don't see any features. There are no error messages.
In an Angular / typescript environment I tried also:
  let polygonStyle = new Style({
  fill: new Fill({
    color: "rgba(255, 0, 0, 0.8)"
  }),
  stroke: new Stroke({
    color: "#ffcc33",
    width: 10
  })
});
let vectorSource = new VectorSource({
  format: new WFS({gmlFormat: new GML3() }),
  url: 'http://osgis.terragis.net/openlayers/examples/gml/polygon.xml',
});
this.vectorLayer = new Vector({
  source: vectorSource,
  style: [polygonStyle]
});

this.map = new Map({
  layers: [
    new TileLayer({
      source: new XYZ({
        url: "https://{a-c}.tile.openstreetmap.org/{z}/{x}/{y}.png"
      })
    }),
    this.vectorLayer
  ],
  view: new View({
    center: transform([0.532597, 47.428810], 'EPSG:4326', 'EPSG:3857'),
    zoom: 12
  }),
  target: "map"
});


texts: How to change the CRS for a Leaflet map at runtime?

texts: I have different layers with different CRS requirements. Since you cannot set CRS per-layer, but only per-map, I need to change the CRS for my map at runtime, depending on which layer is selected. How do I do this in Leaflet?


texts: Identify polygon containing point with R sf package


texts: Changing single character in attribute table with ArcPy?

texts:
I was trying to speed up my data processing, and I was wondering, which would be the python code, that could browse through my shapefiles, and find special single characters (which I define) in attribute table in string fields, and replace them with defined characters.

For example, I have shapefiles, in my "D:\GIS_One" folder. All of the shapefiles have few string(text) fields in their attribute tables. For example, value in some of the fields are  "Day&Night", "Mix!" "Über". In this example, I would like to replace "&" with "_" , "!" with "One" and "Ü" with "U". Just the single characters, not the whole words, since I have another code to change the whole words. I am just looking for a code to do a single character replacement.

I already tried with How to iterate fields and remove Null values and spaces but could not manage to work it out. Also, I know how to to it with Field Calculator, and Excel/LibreOffice, but I am looking for a Python solution so I could speed up my work.

EDIT: As Michael asked, I worked on a first two codes, and now I am working on a third provided from a posted link. Here is working code:

>>> import arcpy
    fc = ["D:\GIS_One\Ride.shp", "D:\GIS_One\Unknown.shp" ]
    fieldList = [f.name for f in arcpy.ListFields(fc) if f.type == "String"]
    if fieldList:
    with arcpy.da.UpdateCursor(fc, [fieldList]) as cursor:
        for row in cursor:
            for i in range (len(fieldList)):
                if not row[i]:
                    row[i] == ""
                elif row[i] == "&":
                    row[i] == "_"
                elif row[i] == "!":
                    row[i] == "1"
                elif row[i] == "Ü":
                    row[i] == "U"
            cursor.updateRow(row)
print "Processing complete"

I get this error:

  Runtime error
  Traceback (most recent call last):
    File "", line 3, in
    File "c:\program files\arcgis\desktop10.1\arcpy\arcpy__init__.py", line 1075, in ListFields
      return gp.listFields(dataset, wild_card, field_type)
    File "c:\program files\arcgis\desktop10.1\arcpy\arcpy\geoprocessing_base.py", line 344, in listFields
      self._gp.ListFields(*gp_fixargs(args, True)))
  IOError: "['D:/GIS_One/Ride.shp', 'D:/GIS_One/Unknown.shp']" does not exist

And if I try with one feature class from database, with this code:

>>> import arcpy
    fc = "D:\GIS_Temp\TEST.gdb\ONE"
    fieldList = [f.name for f in arcpy.ListFields(fc) if f.type == "String"]
    if fieldList:
    with arcpy.da.UpdateCursor(fc, [fieldList]) as cursor:
        for row in cursor:
            for i in range (len(fieldList)):
               if not row[i]:
                    row[i] == ""
                elif row[i] == "&":
                    row[i] == "_"
                elif row[i] == "!":
                    row[i] == "One"
                elif row[i] == "Ü":
                    row[i] == "U"
            cursor.updateRow(row)
print "Processing complete"

I get this error:

  Runtime error
  Traceback (most recent call last):
    File "", line 5, in
  TypeError: 'field_names' must be string or non empty sequence of strings

It would be really helpful if i could manage to make the code to work on shapefiles.


texts: Preparing data for semiovariogram analysis in ArcGIS with large sample sizes?

texts: I am working on a statewide project to create surface water quality grids from point samples.  In some instances, we have large sample sizes (n>50,000) which is problematic when creating semiovariograms in ArcGIS' Geostatistical Analyst.  See my sample locations for groundwater pH.  I'm looking for some literature on how to sub-sample my data prior to interpolation so I can look at the semiovariogram, trends and geostatistics.


texts: open geotiff obtained from NetCDF in python - 'NoneType' object has no attribute 'GetRasterBand'

texts:
I created a geotiff from a netCDF file through the following code, but when I try to access the resulting file in python again, this gives me the following error message: 'NoneType' object has no attribute 'GetRasterBand'  and I do not understand the reason for such! However, the file opens correctly in other software, such as QGIS. Can you help me?

ncdf=netCDF4.Dataset(file,'r')

lat=ncdf.variables["LAT"][:]
lon=ncdf.variables["LON"][:]
var =ncdf.variables["VAR"][:]

(code to convert LAT and LON in X and Y UTM coordinates)

 X=...
 Y=...

 nx=len(X)
 a=numpy.shape(var)
 ny=a[1]

 xres=30
 yres=30

xmin, ymin, xmax, ymax = [X.min(), Y.min(), X.max(), Y.max()]

geotransform = (xmin, xres, 0, ymax, 0, -yres)

 Wkt="""PROJCS["WGS 84 / UTM zone 29N",
            GEOGCS["WGS 84",DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
            PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
            UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],
            AUTHORITY["EPSG","4326"]],
            UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
            PROJECTION["Transverse_Mercator"],
            PARAMETER["latitude_of_origin",0],
            PARAMETER["central_meridian",-9],
            PARAMETER["scale_factor",0.9996],
            PARAMETER["false_easting",500000],
            PARAMETER["false_northing",0],
            AUTHORITY["EPSG","32629"],
            AXIS["Easting",EAST],
            AXIS["Northing",NORTH]]"""

  dst_ds = gdal.GetDriverByName('GTiff').Create('./B1_CATM.tif', ny, nx, 1, gdal.GDT_Float32)
  dst_ds.GetRasterBand(1).WriteArray(refAerosol)

  dst_ds.SetGeoTransform(geotransform)

  srs = osr.SpatialReference()
  srs.ImportFromWkt(Wkt)
  srs.ImportFromEPSG(32629)

  dst_ds.SetProjection(srs.ExportToWkt())

  # AND IN OTHER SCRIPT OPEN THIS GEOTIFF File:

  B=gdal.Open('directoryfile.tif')
  Band=B.GetRasterBand(1)

 # the following error message appears: 'NoneType' object has no attribute
 'GetRasterBand'


texts: Selecting feature programmatically from GeoJSON using OpenLayers?

texts:
I am working on a project that uses OpenLayers (version 2.14)to display a Bing layer(GeoJSON format), I have no problem reading the GeoJSON and display features, but I want to select a feature programmatically, for example, there is a table displaying all the features attributes(GeoJSON format.sample:

  {"type": "FeatureCollection", "features": [
  {"type":"Feature","properties":{},"geometry":{"type":"Point","coordinates":[-7923751.4232522,5233536.7371399]},"crs":{"type":"name","properties":{"name":"urn:ogc:def:crs:OGC:1.3:CRS84"}}}
  ], }

), when I click a row from the table, I want to select or highlight a specific feature on the map using the GEOJSON data in that row.

How can I do that?


texts: Batch convert of viewer styles from QML to SLD

texts:
Until now I was saving styles in QML files.
Now I need them in SLD to be used in GeoServer.
Is there a way how to perform a batch convert?


texts: Virtual fields are lost when duplicate layers

texts:
I've just upgraded to QGIS 2.12.3 - and hoped that the issue which was really time-consuming in 2.12. <3 could have been solved ...

On duplicating layer - with tons of joined tables and virtual fields based on complex expressions - any of these virtually calculated fields are lost.
This issue is already listed on https://hub.qgis.org/issues/14191 - unfortunately just as bug with normal priority.

I think this feature is essential.
are there any workarounds ?

(i know - "Save As" would be a last option. but once more i do not want to duplicate data - which is one of the benefits of join-concept.)

... i have just tried to do "Save Layer as Definition File..."
which seems to work fine concerning the virtual fields
but (!) did "forget" the label expression i've used and one single table join ..


texts: Plot GIS polar data on BlueMarble in Python

texts:
This is my first time working with GIS data and I would like to plot ice sheet extent on a world map with backround centered on the pole.

The shapefile I'm using can be downloaded here: ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/south/monthly/shapefiles/shp_extent/01_Jan/extent_N_197901_polygon_v3.0.zip

I think the coordinates of the data are some sort of polar coordinates because they look like this:

shapefile Reader
    135 shapes (type 'POLYGON')
    135 records (2 fields)
{'type': 'Polygon', 'coordinates': (((-675000.0, 4675000.0), (-625000.0, 4675000.0), (-625000.0, 4625000.0), (-650000.0, 4625000.0), (-650000.0, 4575000.0), (-675000.0, 4575000.0), (-675000.0, 4675000.0)),)}

I'm trying to plot this shapefile on top of a BlueMarble Earth (or any other background). However, the polar coordinates are giving me problems. I tried to plot it using GeoPandas.

import geopandas
import contextily as ctx

df = geopandas.read_file(shp_file)
df = df.to_crs(epsg=3857)

ax = df.plot(figsize=(10, 10), alpha=0.5, edgecolor='k')

This, however, results in several problems as you can see here:

What should I change in this Python script? Or point me in the right direction?

Update:

Based on the comment from Ian Turton I was able to google a bit more and change my code to the following:

df_N = geopandas.read_file(shp_file_N)
df_N = df_N.to_crs(epsg=3031)

fig = plt.figure(figsize=(12,10))
ax1 = fig.add_subplot()

df_N.plot(ax=ax1, figsize=(10, 10), color="white", alpha=0.5, edgecolor='k')
ctx.add_basemap(ax=ax1, url=cartodb_url)
ax1.set_axis_off()

However, the shape and the map are not lining up correctly. If anyone can help me with this step that would be very helpful!


texts: Openstreemap Overpass API to get Couture lines

texts: Is there any way to get contour lines wth OpenStreetMap's Overpass API? I can see them in the cycle map - but I can't seem to get this data at all through Overpass.


texts: Spatially joining small target feature to one of multiple join features using ArcGIS Desktop?

texts:
I am trying to spatially join characteristics from a larger unit (here: county) to a smaller unit (say, local district), such that for every local district I will also know in which county it is in.

The problem arises when a district lies in two (or more) different counties. In this case, I would like the local district to take the attributes from the county that covers most of the district's area. In the screenshot below,

the two local districts marked yellow were assigned attributes from the Greene County despite the fact that most of their area is within the Sumter County. The reason is that I used standard code

target_features = "district.shp"
join_features = "county.shp"
arcpy.SpatialJoin_analysis(target_features, join_features, out_feature_class)

and without any further options specified ArcPy (I guess) simply used the first county information it got hold off.

Is there an option that allows me to account for the area that a district shares with a county?

I have looked at the match_option CLOSEST but I think it does not do the job, because - depending on Green County's total size - its center might very well be closer to a borderline district even if that should be part of Sumter County (like the two yellow districts above).

I am new to ArcPy/GIS.


texts: Using compareRaster (rasters): different extent?

texts:
In R I'm trying to make a RasterStack of RasterLayers, they all have the same format, the only difference is that one of them is a solo raster and the others are bioclimatic variables.
The soil layer had a resolution of ~ 1km and was redesigned to 5km to align with the others.


texts: How do I get the extent of a Gist index?

texts:
I have a bunch of different tables in my database, all of them have a gist index and I want to know of two tables have data in the same area. So what is the easiest way to get the extent of the geometry?

Some tables have a column with points and some tables have a column with polygons..


texts: Grass tools not working in QGIS 2.14.3-ESSEN, OSGEO4W

texts:
Windows 10

QGIS 2.14.3-ESSEN

OSGEO4W install  (current)

64-bit

Processing toolbox.

Grass tools do not work

Neither 6 nor 7. 6 makes sense as I do not see OSGEO4W installing it but 7 tools are not working either.

My Grass folder (7) appears to be C*:\OSGeo4W64\apps\grass\grass-7.0.4*

They are activated (I know they are not in the screenshot).
Grass 7 works as a standalone.

Error message:

  Missing dependency. This algorithm cannot be run :-(  The specified
  GRASS GIS 7 folder does not contain a valid set of GRASS GIS 7
  modules. Please, go to the Processing settings dialog, and check that
  the GRASS GIS 7 folder is correctly configured

What exactly should go in these four yellow boxes?

Now when I try to manually enter my path to grass 7 I get a new error message (below). No msys folder exists in the latest version of OSGEO4W 64-bit so I cannot set this.

.

This solution (QGIS 2.14 standalone GRASS provider missing Msys folder, can't execute tools) appears to point to a spam link now and not mysys.


texts: How to control band order generation in GDAL Merge?

texts:
I am trying to merge a total of 8 GeoTiff files (each one representing a single band) into a single TIF file using Gdal Merge in QGIS. The tool works, but does not allow me to select the order that these layers will be compiled in the output file.

This is problematic as the layer names are NOT conserved in the tool, and therefore it is difficult to tell one band from the other after merging (they are just named Band1, Band2).
I have also tried running this in the CMD (since in QGIS you cannot even try to control layer order), but I couldnt get gdal_merge.py to work. The same for running gdal_merge via PyQGIS.

Here is what I mean by not being able to select layer order in the interface:

The problem is the same across QGIS 2.14 and 2.18.


texts: Optimized walk lists for door-to-door traversal of a region

texts: Given a region in a city, I'd like to make a list of all the street addresses it contains in an order which is reasonably efficient for a person to traverse on foot.  Ideally, I'd like the solution to use free software.  An optimized set of lists for multiple people to traverse simultaneously would be even better.  Does this already exist?


texts: Adding layers to project from a processing tool in QGIS3

texts:
I'm trying to add some layers created during a process to a project.

There seems nothing to be added (TOC empty) but when I close QGIS It says:

  This project includes one or more layers as temporary draft. This
  layers are not saved at the disc and...

Or something like that ( just my own translation)

I defined the layer as:

vl_layerPG = QgsVectorLayer("Polygon", "vl_layerPG", "memory")

Filled with some data and then:

QgsProject.instance().addMapLayer(vl_layerPG, True)

But this is not added to the TOC, almost visible.

Edited 2019/09/16

I tried without "True" first as shown at the pyQGIS cookbook, but It was the same result.

def processAlgorithm(self, parameters, context, feedback):
        dxf_file = self.parameterAsString(parameters,self.INPUT_FILE,context)

        current = 0
        #self.generaCapaPlantaGeneral(datos)

        driverIn = ogr.GetDriverByName('DXF')
        dataSource = driverIn.Open(dxf_file, 0)
        dataSource.ExecuteSQL("SELECT * FROM entities WHERE Layer = 'PG-LP'")
        layerIn = dataSource.GetLayer()
        n_records = layerIn.GetFeatureCount()
        inFeature = layerIn.GetNextFeature()

        geometryPoly = ogr.BuildPolygonFromEdges(ogr.ForceToMultiLineString(inFeature.GetGeometryRef()), dfTolerance = 10)
        geometryPolyWkt = geometryPoly.ExportToWkt()

        ml = QgsVectorLayer("Polygon", "PG_layer", "memory")
        pr = ml.dataProvider()
        pr.addAttributes([
            QgsField("ref",QVariant.String),
        ])

        attrib = []
        attrib.append(46)
        feat = QgsFeature()
        feat.setAttributes(attrib)
        feat.setGeometry(QgsGeometry.fromWkt(geometryPolyWkt))

        pr.addFeatures([feat])
        ml.updateFields()

        QgsProject.instance().addMapLayer(ml)
        ml.updateExtents()

        return {self.MEM_OUTPUT:ml}


texts:
I have a website, http://airport-metars.ddns.net that gets data from FAA and plots on a map using Python and Folium.  The HTML page for the map is created by Python and Folium.  But I would like to add http://vfrmap.com to my map instead.  How can I do that?  I asked this question on openstreetmap.org and they recommended that I ask here.

https://help.openstreetmap.org/questions/71287/adding-overlay-or-tile

I already have an API from vfrmap.com


texts:
I have a map which retrieves Points of Interest from a database and displays them however the Lon Lat need to be transformed into XY to be placed properly.
Currently the all just appear at O,O and I believe they have been transformed to decimal degrees.

I have tried using ol.proj.transform but this returns undefined. Do I have to create a calculation to convert these or is there a method to do this?
The result of this function is returned to a vector layer which loads the features onto a map, albeit in the wrong position.

There are many examples using Openlayers 2 but the syntax is completely different. My non-working example is below:

/* Features conversion */
function wsToCollection(routeData) {

    //console.log(routeData);
    var fc = new Object;
    fc.type = "FeatureCollection";
    fc.features = [];
    for (var i = 0; i < routeData.Points.length; i++) {
        var feature = new Object;
        feature.type = "Feature";
        //console.log(feature);
        feature.geometry = new Object;
        feature.geometry.type = "Point";
        feature.geometry.coordinates = [];
        var Lon = routeData.Points[i].Longitude;
        //var LonTrans = ol.proj.transform(Lon, 'EPSG:3857', 'EPSG:4326'); // not working
        console.log(LonTrans);
        var Lat = routeData.Points[i].Latitude;
        //var LatTrans = Lat.transform('EPSG:3857', 'EPSG:4326'); // also not working
        feature.geometry.coordinates[0] = routeData.Points[i].Longitude;
        feature.geometry.coordinates[1] = routeData.Points[i].Latitude;
        feature.properties = new Object;
        feature.properties.name = routeData.Points[i].Description;
        feature.properties.Guid = routeData.Points[i].Guid;
        feature.properties.LongId = routeData.Points[i].LongId;
        fc.features.push(feature);
    }

    return fc;
}


texts: Calculating distance of point from arc/polyline using PostGIS?

texts:
Using PostgreSQL/PostGIS, I am trying to query the database for the closest arc to a particular GPS point. Each Arc is made up of two nodes that are GPS points on a map. Is there a way I could query this using Postgres and also get the distance of the point from the arc?

I have been trying to find the road to which a GPS point is closest to. I have a postGIS database which i loaded with the map of the road network.
According to a previous question asked here Calculate distance of point from an arc/polyline this was to be done be creating a straight line then look for the shortest distance form this line. The road network is however not straight thus this wont work.

Anybody with a solution?


texts: Globbing in Windows to process images - how do bash loops look in DOS?

texts:
I'm doing some batch image processing with imagemagick on OSX using this line:

for file in $(cat all.txt); do convert -crop 50X50@ +repage +adjoin "$file" "$file"_tile_50x50@_%d.jpeg; done

It takes a filelist as an input, chops each file into 50x50 pixel tiles, and labels the output with the original filename and an index number (%d). The quantity of imagery is about 400+ source images, and hence 1,000,000 output images.

I have to move the processing to a Windows 10 box. Can someone help me translate the above loop across, or suggest a better way to achieve the same outcome?


texts:
I need to convert the 1/3 arc second DEM files from the NED into WGS84 ellipsoid heights. The documentation says that the NED is referenced to NAVD88. I've tried using gdalwarp as follows:

gdalwarp imgn45w124_13.img imgn45w124_13_wgs84.tif -t_srs "+proj=longlat +ellps=WGS84"

But this doesn't change the heights at all, and they are still off. What is the correct way to do this?

I've also tried the following without success just incase GDAL isn't able to read in the vertical datum info.

gdalwarp imgn45w124_13.img imgn45w124_13_wgs84.tif -t_srs "+proj=longlat +datum=WGS84" -s_srs "+proj=longlat +datum=NAVD88"


texts: WMS: How can I connect to deegree (localhost)?

texts:
I have deegree version 3.4.16 running on 127.0.0.1:8080 (localhost). I want to connect QGIS version 3.10.5-A Coruña to deegree (localhost).
I have set the following parameters in "Create a New WMS/WMTS Connection":
Name = deegree2
URL = 127.0.0.1:8080 (or localhost:8080, both tried)
The following error message appears:
"Failed to download capabilities. Download of capabilities failed. Protocol "" is unknown
What is missing? What went wrong


texts: Selecting features with field not contained in all records of another field using QGIS?

texts:
I am quite new to QGIS 3.4.4 and it is complicated at times to get the proper thinking style.

I want to shift the location of some points in my vector layer.
I have a vector layer with points of waterbodies ("waterpoints") with a unique identifier field OBJECTID. I have another vector layer with some points which I relocated on the field. They have a field "Name" which refers to the OBJECTID.

I thought if I can delete all the original points whose "OBJECTID" is contained in "Name" of the other layer and thereafter merge the layers I would be fine. But I can´t make a selection or an export which can check that.

First I didn't manage to get the field from another layer, so I just merged them to run the expression within one layer.

I tried:

OBJECTID IN  (Name)

but it does not work...

I really don´t know how to proceed.
I would like to know how to properly write this expression or if there is an alternative method for relocating those points.


texts: Test for no selection without a definition query

texts:
I am trying to put together a little python script to check for whether there is a selection. The dataset I have is not robust enough to have an definition query on it, so I want it to just check whether there is a selection. arcpy.GetCount_management does this, and when I use the Describe function it does this too, but whenever I don't have a selection it is giving me the total number of rows in the dataset. From what I've read in the literature, the Describe function at least should return a zero value, because none are selected. Can anyone tell me why I can never get a zero value from either of these tools?

This is the piece of code I've been playing with most recently

len(arcpy.Describe("Footpath_centreline_Layer").FIDSet.split(';'))

so when I have nothing selected, 0 out of 16 records - I get a return of 16, and I need for it to be 0. I'm not sure how to achieve this, given my very limited python knowledge.


texts: Running into Java Heap Size issues with OpenTripPlanner

texts:
I'm trying to load up OpenTripPlanner on my PC and utilize it.

I get everything set to run, then it gets about two minutes into the OSM file nad icks out "exception in thread "main" java.lang.OutOfMemoryError: Java heap space", followed by "at", then listing a whole mess of places it's out of memory.

I've tried -Xmxing every file size it will allow, to no avail.

(I'm running Win8.1 with 16gb RAM).

I don't have access to a shell running Java to use this with, so hopefully I'll be able to make my system work.


texts: I have a single band raster image with values of pixel 0 and 1. There's another line vector layer over the same area. I need to find the percentage of line which is overlapped by the pixels of value 1. I am thinking of using clipper in qgis or mask function in R but can't exactly figure out the way. Is there any way of doing this in QGIS or R?


texts: Calculate 70% of the area inside a polygon based on the polygon centroid

texts:
I am trying to generate an area within a polygon based on the polygons centroid that covers 70% of the containing polygons area. I generated the polygon centroids as points but I do not know how to generate a buffer based on both the polygon centroid as the start of the buffer location and constraining the desired polygon to be 70% of the original polygon size.

FYI - polygons are parcels


texts: How much process memory is needed for ogr2ogr to process a GeoJSON file? -

texts: I am trying to evaluate how many resources ogr2ogr needs to process a GeoJSON file; there's a concern that it's using much more memory than we would expect.  All we are doing is converting a GeoJSON file, filled with point features, to create another GeoJSON file - to validate the file, and convert the coordinate system to our target CRS.   When I run it under Ubuntu Linux, and run the TOP command, the process seems to take about 134000 KB  of virtual memory, regardless of the file size.  However, the amount of buff/cache used seems to vary with the size of the file being converted.  How much memory should a process have to run ogr2ogr?


texts: How to add GeoTagged images to the map and store as KML

texts: I have the requirement that i need to add a GeoTagged Image into the map and save that image in KML format. Give the approach to do this task.


texts: Aggregating administrative boundaries using ArcMap

texts:
I am working with a shapefile with polygons at the state level which I am grouping using the dissolve tool.

How can I save the resulting layers produced from the several dissolve output layers into one layer that retains the new aggregated properties of the regions?

Merging does not retain the new boundaries.


texts: PostGIS: How to ST_Clip() a raster with a polygon?

texts:
I have done:

Import a SHP file into PostGIS by PostGIS SHP loader GUI.
Import a raster into PostGIS:

\raster2pgsql -I -s 4326 landsat.TIF -t 500x500 public.landsat > query.sql

After importing all into postGIS, I have tried exporting raster to PNG by ST_AsPNG. It run well. Besides, I run the example of ST_Clip() in PostGIS docs. It worked.

Question: I am trying to clip a raster with a polygon of the SHP file. It returns nothing.
I tried ST_Intersects(rast, geom polygon) like this. It failed.

Can anyone help me on clipping a raster with a pylygon of SHP file?

Viewing in QGIS, the raster overlaps the polygon. I'm using PostGIS 2.0.

EDIT: Here's what I do with ST_Intersects()

SELECT
ST_Intersects(polygon.geom, raster.rast)
FROM
  public.vnm_adm2 as polygon,
  public.landsat as raster;

And it returns false in all rows.

UPDATE:  There're some things to be noted here:

I used ST_AsPNG(), which remove spatial reference to output file.
Switch -s [4326] doesn't work yet (PostGIS 2.0), so the rasters ain't be converted to polygon's EPSG, then it returns false in all rows.


texts: How to align bounding box when CRS match

texts:
I have two layers, US State boundaries and international boundaries for all countries.  I know that due to different data sources and resolutions the lines won't match exactly but I'm getting a strange offset on the two layers (see below)...

the CRS match but the bounding boxes are different...  How do I align these?

US States shapefile download

Dept of State Large Scale International Boundaries

# Load packages
library(dplyr)
library(sp)
library(rgdal)
library(raster)
library(rgeos)
library(ggplot2)
library(ggmap)

# Read in and view GIS layers
# Vectors
lsib1 <- readOGR(dsn = "R:\\myGIS\\DoS\\Department_of_State_Large-Scale_International_Boundary_28LSIB29",
                 layer = "Department of State Large-Scale International Boundary (LSIB)")
states1 <- readOGR(dsn = "R:\\myGIS\\Census\\tl_2017_us_state",
                   layer = "tl_2017_us_state")

# Project layers to Albers
crs(states1)
crs(lsib1)
lsib2 <- spTransform(lsib1, proj4string(states1))
crs(lsib2)

# Visualization
par(mfrow = c(1,1))
plot(states1)
plot(lsib2, add = TRUE)

In response to comment: these layers align just fine (with minor incongruities) in GIS:

Above is my complete code.


texts: SQL filter to show one layer only if another layer is missing

texts:
I have two oracle tables that represent buildings in my town. The first table represent buildings that are surveyed and measured accurately but doesn't cover all the parts of the town while the second table represent not so accurately measured buildings but covers all the town.

table_1 - accurate data - partial coverage
table_2 - not accurate data - full coverage

I want to write some sort of an SQL filter to show the data from table_2 only if no data in table_1 are present in that area. In other words, I want the whole town to be covered by buildings and to show all the areas that have the accurate buildings and fill up the rest of the areas with the not so accurate buildings!

The tables are stored in a Oracle 11g database. The map application I am using is built on MapServer and the layers are basically defined in XML files that accepts SQL filters.

So far, I have done the following:

select table_2.* from table_2, table_1
WHERE SDO_ANYINTERACT(table_2.geom, table_1.geom)= 'TRUE';

The problem with this method is that it gives me the places of the inaccurate buildings when it interacts with the accurate buildings. What I really want is to set the SDO_ANYINTERACT to 'FALSE' to get the inaccurate buildings in the places where no interaction is occurring between the buildings but of course Oracle gives error when setting that to 'FALSE'.

Any suggestions?


texts: Referencing fields in code block of Field Calculator in ArcGIS for Desktop?

texts: I'm trying to use values from a field to fill a new field with a truncated version. The 'str.rindex' function is working out well, but I can't seem to figure out how I'm supposed to reference the field I need truncated. I found a few answers on here explaining that the !field! convention won't work, but no solutions for a workaround.


texts: How can I customize GLMap offline vector map?

texts:
I am using GLMap from https://getyourmap.com/ and I have two questions:

In my vector map the street name is shown in two languages I want to show in one language. How can I change that?
Where can I give information about adding marker or set locations?


texts:
I am unable to open any new map. Add Data doesn't show any basemaps or Arcgis Online Maps. As soon as I try to add any of the base maps, my program closes on me.

Yet I am able to work on all my old saved maps with basemaps in them so far. Also, I can work with my database without adding any basemap or arcgis online maps.

To make a new map, the round about way I have figured out is to add a basemap from ArcCatalog/GIS servers/arcgis on services.arcgisonline.com

But I want my Add Data to add basemaps to my new maps.

How do I do that?


texts:
Does anyone know if there is a way to track or get email notifications if someone has made changes to attributes via online mapping?
Our services are hosted through Arcserver 10.1 and an SDE database.

Currently we are using Geocortex Essentials for our online mapping needs, however I'm wondering if there this capability exists across the board?

I would just like some way of tracking changes and data management.

Thanks.


texts: Update/refresh edited PostGIS table in QGIS

texts:
I imported a PostGIS table (with geometry) from my PostgreSQL database as a layer in QGIS.
Later, I edited the table using PostgreSQL UPDATE commands.

Is there a way to quickly update (like a refresh ↻ command?) the layer in QGIS with the new attributes or geometry without having to remove the layer and reconnect with the database?


texts: Trying to create a new coded domain and Arcpy is rejecting my domain name

texts:
For some reason I'm getting a "ERROR 000192: Invalid value for Domain Name" in the below script. The domain name is "Road_Class_Types". I can't figure out why it's giving me this error. Why doesn't it like it? What am I missing?
# Import modules
import arcpy
from arcpy import env

# Set the workspace (to avoid having to type in the full path to the data every time). Update foe each new file.
arcpy.env.workspace = r"C:\Program Files (x86)\ArcGIS\Desktop10.8\ROBERT FILES\EXTRACTION WORKING FOLDER"

# Set parameters
domName = "Road_Class_Types"
gdb = "ROBERT_WORKING_EXTRACTION.gdb"
inFeatures = "ROBERT_WORKING_EXTRACTION.gdb\TransportationGroundCrv_AddFieldTest"
inField = "AP030_Road_Class"
fieldName1 = "AP030_Road_Class"
fieldPrecision = 9
fieldAlias = "Road_Class"
fieldLength = 200

# Process: Create the coded value domain
arcpy.CreateDomain_management("ROBERT_WORKING_EXTRACTION.gdb", "Road_Class_Types", "Road Class types per LTDS schema", "TEXT", "CODED")

# Store all the domain values in a dictionary with the domain code as the "key" and the domain description as the "value" (domDict[code])
domDict = {"Class 5 Road": "Class 5 Road - Unimproved (ZI016_ROC = 1), Fair Weather (ZI016_WTC = 2), one (1) or two lanes (2) or less, (LTN = 1 or 2, FCO =7", \
           "Class 4 Road": "Class 4 Road - Flexible Pavement (ZI016_ROC = 3), Limited All-weather (ZI016_WTC = 4), one (1) or two lanes (2) or less, (LTN = 2), FCO =7", \
           "Class 1 Highway": "Class 1 Highway - Rigid Pavement (ZI016_ROC = 7), Divided Highway, All Weather (ZI016_WTC = 1), three lanes (3), (LTN = 3), FCO = 5 or 6", \
           "Class 2 Highway": "Class 2 Highway - Rigid Pavement (ZI016_ROC = 7), Divided Highway, All Weather (ZI016_WTC = 1), two lanes (2), (LTN = 2), FCO = 5 or 6", \
           "Class 3 Highway": "Class 3 Highway - Rigid Pavement (ZI016_ROC = 7), Undivided Highway, All Weather (ZI016_WTC = 1), one (1) or two (2) or less, (LTN = 2), FCO = 7", \
           "Superhighway": "Superhighway - Rigid Pavement (ZI016_ROC = 7), Divided Highway, All Weather (ZI016_WTC = 1), four or more lanes (4), (LTN>=4),FCO = 5 or 6", \
           "noInformation" : "noInformation"}

# Process: Add valid types to the domain
#use a for loop to cycle through all the domain codes in the dictionary
for code in domDict:
arcpy.AddCodedValueToDomain_management(gdb, domName, code, domDict[code])

# Process: Constrain the material value of distribution mains
arcpy.AssignDomainToField_management(inFeatures, inField, domName)


texts: Processing framework input variables from QGIS-2 not working in QGIS-3

texts:
A simple script created within the QGIS 2.18 processing toolbox:

##example=string
print example

Upon running would produce and input box in a window with a run button:

Once executed would produce the expected print in the python console

When testing the same code in QGIS-3.0, changing the print to python 3's print(), I get an error where the input section is seemingly ignored as a comment:

Does the processing framework no longer accept inputs in this way?


texts: how to preserve color-interpretation when merging TIFFs and adding alpha band

texts:
We have a collection of black-and-white map tiles as *.tif images with accompanying *.tfw files, that we'd like to merge into one single large map and cut with a polygon (turning the area outside the polygon transparent).

gdalinfo output for them looks like this:

Driver: GTiff/GeoTIFF
Files: 702_232.tif
       702_232.tfw
Size is 4000, 4000
Coordinate System is `'
Origin = (702000.000000000000000,233000.000000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Metadata:
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=96
  TIFFTAG_YRESOLUTION=96
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
  MINISWHITE=YES
Corner Coordinates:
Upper Left  (  702000.000,  233000.000)
Lower Left  (  702000.000,  232000.000)
Upper Right (  703000.000,  233000.000)
Lower Right (  703000.000,  232000.000)
Center      (  702500.000,  232500.000)
Band 1 Block=4000x8 Type=Byte, ColorInterp=Palette
  Image Structure Metadata:
    NBITS=1
  Color Table (RGB with 2 entries)
    0: 255,255,255,255
    1: 0,0,0,255

So this seem to be 1-band TIFFs with a color palette.

To merge them, we're first building a VRT with

gdalbuildvrt -a_srs 'EPSG:21781' merged.vrt *.tif

(as you can see, we also set the source SRS, which is missing in the input data)

The individual input files (the tiles) and the VRT display fine (black and white, as intended) in QGIS 2.18.

To build an actual single image file (and at the same time cut the result to the polygon area1), we then use

gdalwarp \
    -cutline swissBOUNDARIES3D/FGDB/swissBOUNDARIES3D_1_3_LV03_LN02.gdb \
    -cl TLM_KANTONSGEBIET \
    -cwhere "NAME='St. Gallen'" \
    -dstalpha \
    \
    -co TILED=YES \
    -co COMPRESS=LZW \
    -co BIGTIFF=YES \
    \
    --config GDAL_CACHEMAX=40% \
    -multi \
    -et 0 \
    -wo NUM_THREADS=ALL_CPUS \
    -wm 5000 \
    \
    merged.vrt  merged_and_cut.tif

The result is a 2-band TIFF (the 2nd band being alpha, required for the transparency after cutting) without a color palette. When displayed in QGIS, it appears in bright green instead of black and white.

If we add -setci to the options for gdalwarp we get a 2-band TIFF with a (sensible looking, even though much larger) color palette. But in QGIS 2.18, it also appears in bright green by default, at least until go to layer properties → style and manually set "Render type" to "Palette".

As the original files and the VRT work fine, I wondered whether we can use the  "Raw mode" documented at http://gdal.org/frmt_gtiff.html to preserve the first band's color interpretation, but there's no indication on how to use with gdalwarp. (Replacing merged.vrt by GTIFF_RAW:merged.vrt in the the gdalwarp invocation results in an error about the file not being found.)

How can I preserve the color interpretation with gdalwarp? (So that it is correctly understood by QGIS and other software.)

Example data

For licensing reasons I think I cannot share the tiles I want to merge, but I've produced these two example tiles (with nonsensical content) that have the same properties and show the same behavior as the original tiles:

Inputs (TIFFs and TFWs): ownCloud directory | download as ZIP archive
gdalbuildvrt -a_srs 'EPSG:21781' merged.vrt ../input/*.tif results in merged.vrt (2 KB) (direct download)
gdalwarp -dstalpha -co TILED=YES -co COMPRESS=LZW merged.vrt merged.tif results in merged.tif (17.6 MB) (direct download)

1 Thus why we use gdalwarp. If we wouldn't cut to a polygon, gdaltranslate instead of gdalwarp would suffice to materialize the VRT to a single TIFF.


texts: QGIS 2.0.1 linking table issue, MSSQL 2008 database

texts: I'm trying to link tables in QGIS 2.0.1 and MSSQL 2008 database.  I can link them but when I try to open the attribute table QGIS throws an error and shuts down creating a *.dmp file.  I saw that table joining when using databases that QGIS does not support MSSQL but the funny think is that QGIS 1.8 allows the join and attribute tables to work.  Anyone have any information regarding this?  I haven't been able to find anyone else having this problem online.


texts:
I'm using Visual Studio 2010 to open and build the GDAL 1.10 source (using the "makegdal10.sln"). This builds without an issue.

Now I want to include the ECW support by including the ECW v5.0 SDK (downloaded from Intergraph). It is installed to "C:\Intergraph\ECWSDK" and the relevant nmake.opt now looks as follows:

ECWDIR  =  "C:\Intergraph\ECWSDK"
ECWFLAGS = -DECWSDK_VERSION=50 \
       -I$(ECWDIR)\include \
       -I$(ECWDIR)\include/NCSECW/api -I$(ECWDIR)\include/NCSECW/jp2 \
       -I$(ECWDIR)\include/NCSECW/ecw
ECWLIB  =  $(ECWDIR)\lib\vc90\win32\NCSEcw.lib

However when building I get the following errors and have no idea what it means:

Error   6   error LNK1120: 2 unresolved externals   (file)gdal110.dll

Error   5   error LNK2019: unresolved external symbol "__declspec(dllimport) public: __thiscall std::basic_string<wchar_t,struct std::char_traits<wchar_t>,class std::allocator<wchar_t> >::~basic_string<wchar_t,struct std::char_traits<wchar_t>,class std::allocator<wchar_t> >(void)" (__imp_??1?$basic_string@_WU?$char_traits@_W@std@@V?$allocator@_W@2@@std@@QAE@XZ) referenced in function "private: void __thiscall ECWDataset::WriteHeader(void)" (?WriteHeader@ECWDataset@@AAEXXZ)   (file)ecwdataset.obj

Error   4   error LNK2019: unresolved external symbol "__declspec(dllimport) public: wchar_t const * __thiscall std::basic_string<wchar_t,struct std::char_traits<wchar_t>,class std::allocator<wchar_t> >::c_str(void)const " (__imp_?c_str@?$basic_string@_WU?$char_traits@_W@std@@V?$allocator@_W@2@@std@@QBEPB_WXZ) referenced in function "public: virtual enum CPLErr __thiscall ECWDataset::SetMetadataItem(char const *,char const *,char const *)" (?SetMetadataItem@ECWDataset@@UAE?AW4CPLErr@@PBD00@Z)   (file)gdal-1.10.1\ecwdataset.obj

Error   8   error MSB3073: The command "cd C:\Will\Mapserver\gdal-1.10.1\ && nmake -f makefile.vc MSVC_VER=1500 DEBUG=1 && nmake -f makefile.vc MSVC_VER=1500 DEBUG=1 install" exited with code 2.  C:\Program Files (x86)\MSBuild\Microsoft.Cpp\v4.0\Microsoft.MakeFile.Target

Error   7   error U1077: '"C:\Program Files (x86)\Microsoft Visual Studio 10.0\VC\bin\link.EXE"' : return code '0x460'  ...\gdal-1.10.1\NMAKE   makegdal


texts:
I'm currently working to find the intersections between linears in a shapefile. I'm a newbie to the GDAL/OGR library. I've done some googling and came up with what seems to be a useful example(see below). There is a problem though. I don't see the method ExportToWkt(). My intellisense(or whatever it's called in IDLE) for the variable "intersection"doesn't have that option. Can anybody tell me what I'm doing wrong. Or just give me some general guidance to comparing two features in a shapefile and finding the intersection(if they intersect at all). Thanks.

from osgeo import ogr

wkt1 = "POLYGON ((1208064.271243039 624154.6783778917, 1208064.271243039 601260.9785661874, 1231345.9998651114 601260.9785661874, 1231345.9998651114 624154.6783778917, 1208064.271243039 624154.6783778917))"
wkt2 = "POLYGON ((1199915.6662253144 633079.3410163528, 1199915.6662253144 614453.958118695, 1219317.1067437078 614453.958118695, 1219317.1067437078 633079.3410163528, 1199915.6662253144 633079.3410163528)))"

poly1 = ogr.CreateGeometryFromWkt(wkt1)
poly2 = ogr.CreateGeometryFromWkt(wkt2)

intersection = poly1.Intersection(poly2)

print intersection.ExportToWkt()


texts:
I am working with raster/brick data in R and running into memory issues. I am basically trying to use raster objects to perform some simple agent-based modeling, but the objects keep increasing in size for reasons I can't understand.

Below is some code that demonstrates what I am talking about.

#### Test memory ###
library(pryr)
library(raster)

br = brick()
for(i in 1:25){
  r1 <- raster(nrow = 100,ncol = 100)
  r1[] <-rbinom(length(r1),1,.5)
  br = addLayer(br,r1)
}
out = NULL
for (i in 1:nlayers(br)) {
  br[[i]][] = 0
  print(object_size(br))
  out[i] = object_size(br)
}

out = cbind(1:25,out)
plot(out[,1],out[,2])
lm(out[,2] ~ out[,1])

To me (limited understanding of R memory usage), the object size shouldn't change. I am simply converting 1s to 0s.

Can anyone help me understand this and find a workaround?


texts: Creating numerical mask file with elevation points of area with QGIS?

texts:
I'm using QGIS and I have to create a mask file (.txt) as an input for a hydrological model. The file should be a numerical grid map that represents the elevation points of a catchment area, like the example in the following picture (where values 0 represent points which are out of the interest zone).

I have a 40m DEM of the region and the shapefile of the catchment area. I need a cell size of 100m and the points out of the area must be -9999.

Does anyone know how I can do this with QGIS?


texts: QGIS Contour plugin: Creating complex isolines leads to isolated areas

texts:
I use the contour plugin (v1.1) within QGIS 2.4 to create contour lines. I have streets with different travel costs. I observe problems for streets with low travel costs if they are surrounded by streets with high travel costs. The low cost street is "piercing" into the low speed area (cf. picture). The points are created as the travel cost markers from Dijkstra algorithm using pgrouting 2.0.5. The Dijkstra algo. calculates correct values for every street.

The problem I observe is now, that due to low speed areas left and right to the highway, the isolines are connecting both areas over the highway (yellow). Based on values, the correct behavior would be the one shown in green.

Can someone recommend me a way of preventing such results?

Before I used GDAL TIN interpolation and contour mapping, which leads to better results if it comes to complex networks. On the other hand it adds transformation errors (rasterizing) and extra steps to the calculations.


texts: Publish geoprocessing service on ArcGIS for Server?

texts:
I am following the instructions described here and I have some enquiries about publishing a GP service.

First of all, how should I publish a GP service on a production environment where the database is clean, without any data. Should I use random input parameters so as to produce a result (which will be empty) , and in turn select Share As > GP Service ?
I cannot use SD file for deploying, So I have to do it using ArcMap


texts: Is there an EPSG code appropriate for locally projected data?

texts:
I often use GIS to track some very local things. Particularly, mapping on microscope slides works quite well. These obviously cannot be georeferenced in any meaningful way.

Perhaps the right thing to do is just not to set a CRS, but I'd like to remove ambiguity any way I can. Any ideas?


texts: vector tool missing when install QGIS 1.7.4 on new computer

texts:
I installed QGIS version 1.7.4 a few months ago. Under the menu Vector > Data Management Tools I have two options to make a join: 'Join Attributes' and 'Join attributes by location' (see screenshot).

My colleague installed the same version a couple of weeks ago.  In here version, the option 'Join attributes' is missing.  Today I installed QGIS 1.7.4 on some computers in a classroom and the same option was missing.

It is a handy tool to join a shapefile with a DBF, instead of doing it through Layer properties > tab Joins.  Can I find it somewhere else?


texts: Doubts about how pgRouting topology is built and how flexible it is

texts:
I'm learning pgRouting, to validate if it can be used in a project. I have several doubts. Mostly related with how could I "fill the gap" between the "bricks" it provides and the real world applications based on routing I normally use. I will summarize the questions in a few ones. Sorry if they're too silly. I have little experience with topology.

EDIT: PolyGeo wisely suggests that I focus this question in one important topic, and leave the rest for different questions. I will reformulate my doubts to create a general question important for me. The already existent responses will still make sense. Anyway, I keep the original questions at the end.

My general doubt is about the topology creation process and how flexible is.

If my data consists in a bunch of lines, and I tell pgRouting to create a topology from that, is a node created at any intersection, apart from the nodes created at the beginning and end of any line? And what about intersections at different level? (for example: a bridge above a motorway: http://goo.gl/pVeE4j, that's not really an intersection) How are these handled? --> EDIT: I've already found this http://goo.gl/ojGmuH. I think it is a response. I need to investigate it.

Once the topology has been built, what do I do if I need new nodes for a query? Example: I want the shortest path from any arbitrary point in my route to another arbitrary point, not just from one existent node to another one (example: the Google Maps feature "give me directions from my current location to this marker I've just put in the map")

Regarding the topology, my last doubt is if the PostGIS topology and the pgRouting topology are related in any sense.

These are the original 6 questions I made:

About creating the topology: you feed the engine with a few lines, as the base of the topology. I guess a new vertex is created every time two of these lines intersect, apart from the vertex created at the beginning and the end of each line. Am I right?
When two lines intersect at different level, how is it handled? (if it's really handled). For example, a bridge above a motorway. Something like this: http://goo.gl/pVeE4j
I've seen several examples using really simple data: a few lines forming a simple network. I've also seen examples using OSM data imported with osm2po or osm2pgrouting. OSM data is great, and complex. But what if I want an intermediate point? Not as simple as a single table with a few dozens of prepared linestrings and not as complex as OSM datasets. Sometimes, we need to work with data provided by public institutions (old digitized maps, cadastre datasets, etc). Is there any way to integrate data from different sources? Any method to ease the compliance of my data with the minimum requisites for pgRouting to create a topology based on it? So, to go from chaos (different sources, different methods to obtain the data...) to a polished set of PostgreSQL tables, ready to work with?
The minimum path algorithms always need two nodes (from A to B). No problem if those points are actual vertices of my topology. But, what if I want to calculate a path from this exact point where I am right now to another arbitrary point? Like the Google Maps calculation: "give me the path from my current location to X". Creating new vertices every time an arbitrary path is requested and rebuilding the topology sounds like nonsense for me.
About the turn restrictions, I find the underlying logic difficult to understand. What do the parameters of restriction SQL mean? In the example, a restrictions table is created. Why those fields?
Do pgRouting topology and PostGIS topology know each other? If not, are they planning to converge/collaborate in any point?

So, summarizing, in my opinion pgRouting looks like a really great framework, but I think there are a lot of things to learn, understand and build in order to construct a routing engine using it.


texts:
I installed GRASS GIS 6.4.3 on Ubuntu Precise 12.04 by sudo apt-get install grass grass-doc, as recommended by the Wiki. It runs with the wxpython GUI, but when I use any command (like r.in.gdal or v.digit etc.) from the command line, i get the exception:

Traceback (most recent call last):
    File "/usr/lib/grass64/etc/wxpython/gui_core/forms.py", line 69, in <module>
    from core import globalvar
    ImportError: cannot import name globalvar

I opened the python interpreter (actually running python 2.7.3) and tried to import core; that works but, in fact, there's no module globalvar under core.
I am able to use all grass commands from the GUI, but the same commands raise the ImportError exception if run by the command line.
This happens only when starting GRASS with the wxpython GUI, not with the Tcl/Tk one. I assume that GRASS looks for globalvar in the main core package rather than in the wxpython.core one.


texts:
I have had variable success with photo2shape in QGIS. I can't figure out what is the success/no success prerequisites. Any thoughts?
Currently I have a directory with all 334 photos having geotag data. QGIS stalls when photo2shape is set to work, and QGIS is hard to kill.
On a test of 5 of those photos it runs fine.