texts:
I am an architecture student with no experience with QGIS. However, I need some .dwg drawings of a town in Italy called Certaldo. I found detailed drawings for the town and had to download it as a .qgs file. I have the file open in QGIS and from the little reading I have done, the layers I am trying to export needs to be a vector layer in order to export it as a .dwg file. However it seems to be a raster layer.

How do I convert this raster layer into a vector layer?


texts: How to set version of wfs-t in openlayers-3?

texts:
I'm using wfs-t in OpenLayers-3.6 according to here. I modified a feature and then sent a wfs request as follows:

var properties = modifiedFeature.getProperties();
delete properties.boundedBy;
newFeature = new ol.Feature(properties);
newFeature.setId(features[i].getId());

var node = wfs.writeTransaction(null, [newFeature], null, {
    featureType: "Calk:Point", //geoserver Layer Name
    featureNS: "http://itsme.calk.ir", // Edit Workspace Namespace URI
    srsName: "EPSG:4326"
});
var s = new XMLSerializer();

var url = "http://localhost:8084/geoserver/wfs";

Ext.Ajax.request({
    url: 'proxy.cgi?url=' + encodeURIComponent(url),
    method: 'POST',
    xmlData: s.serializeToString(node),
    headers: {
        'Content-Type': 'text/xml'
    },
    success: function(response_) {
        response = response_;
    }
});

The default version that OpenLayers works is "1.1.0", In this setting, I want to work with version "2.0.0". How can I set the version?


texts: Calculate geographic size in meters given lower left coordinate and upper right coordinate

texts: How do you find the geographic size in meters given a lower left coordinate and upper right coordinate. The projection would be WGS84 projection. Ex. 33,-119,35,-117. What would be the area in meters that covers?


texts: Querying highways with Overpass QL for calculating routes

texts:
I want to query highways with Overpass to calculate routes with Graphhopper. With some trial and error I came to the following query and it seems to work (tested in overpass turbo):

(
  ._;
  way[highway]({{bbox}});
  node(w);
);
out;

I still have trouble understanding what exactly the query does, especially the default set ._;. Is it actually redundant here?

Also, I really only need to calculate routes for pedestrians and cyclists. What would be the best approach there? Expanding the query above with filters or just limit the vehicle in Graphhopper?


texts: How can I make an OSM map without icons and labels?

texts:
I want to create a web application that include Open Street Map .
I use OpenLayers library to create map , but now I want to create map without icons like hospitals,coffeshop and other icons.

This images explain:


texts: Creating schematic view of line feature in ArcGIS Desktop?

texts:
I have a line feature (polyline with curves) in ArcGIS Desktop, and I want to create a linear representation (straight line) of the feature, furthermore I want to be able to locate events in the linear representation. This is something like I would like:

Making some research I found something similar called "column view" unlucky my PC rebooted an does not saved the history of the browser... so I lost the document. (I have been searching for it again, but I have not been able to find it)

In the petroleum industry they call it "Alignment Sheet", is a representation of the objects that are located over a route but in a straight line. But keep in mind that is not linear referencing. Think in this way: you have a road, with its events, it is a curve, the idea is to straight it and view the event over that line.


texts:
Where can I download the exercise data for The Free Quantum GIS Training Manual?
The pdf of the manual just mentions the sources of the data, but doesn't give any links.


texts: How to add labels to a WMS layer using OpenLayers and GeoServer?

texts: I have point of data in mysql server. I am displaying this layer as WMS using OpenLayers and GeoServer. How can I add label to this layer?


texts: I have GIS points (x, y data) for every minute over a landscape. If I load them up as a CSV into QGIS they all show up as points, which you would expect. However, what I would like to do is use these points to create a grid. In other words, each point is a corner of a 1x1 minute grid. How can I do this?


texts: Is it possible to route shapefiles using python and without ArcGIS, QGIS, or Pgrouting?

texts:
The question is pretty much in the title.

Is there a way to route through a road line shapefile using python and some library without having to resort to converting the data into PostGIS/pgrouting or going through desktop software API such as ArcGIS and QGIS?


texts: Making loop over shapefile attribute using ArcPy?

texts:
I have a polyline shapefile. The shapefile has a field called "NAME"

I would like to make a loop over each line (A,B,C,D,...).

This is the python code for only first line named "A"

profileLine = WORK_DIRECTORY + r'\CrossSections_onlyA.shp'
profileTargets = WORK_DIRECTORY + r'\DTM.tif'
OutTable = WORK_DIRECTORY + r'\OutTableA.dbf'
arcpy.StackProfile_3d(profileLine, profileTargets, OutTable)

I don't know how to do it over all of them with a loop. Also I would like to save the table with the name of each Line. for example: OutTableA.dbf, OutTableB.dbf, OutTableC.dbf and so on.


texts: Unable to create tirex gdal datasource for OpenTopoMap

texts:
I try to set up an OSM tile server using tirex on a server running Ubuntu server 16.04 LTS 64-bit.
In the style definition (opentopomap.xml) there are several gdal data sources definied. E.g. the following color relief layer:

<Layer name="relief-5000">
<StyleName>relief-5000</StyleName>
<Datasource>
    <Parameter name="type">gdal</Parameter>
    <Parameter name="file">dem/relief-5000.tif</Parameter>
    <Parameter name="format">tiff</Parameter>
</Datasource>

At the relative path to the style-file (dem/relief-5000.tif) I have the tiff file with the following access rights:

-rwxr-xr-x 1 postgres postgres 8866 Feb  5 01:29 relief-5000.tif*

My tirex user is in the postgres group and should be able to read the layer's file

uid=112(tirex) gid=119(tirex) groups=119(tirex),118(postgres)

But when I start tirex and check /var/log/daemon.log I get the following error message:

Feb  5 10:27:23 ubuntu-server tirex-backend-mapnik[4739]:
Renderer started (name=mapnik)

Feb  5 10:27:23 ubuntu-server tirex-backend-mapnik[4739]:
cannot add /etc/tirex/renderer/mapnik/opentopomap.conf

Feb  5 10:27:23 ubuntu-server tirex-backend-mapnik[4739]:
Could not create datasource for type: 'gdal'
(no datasource plugin directories have been successfully registered)
encountered during parsing of layer 'relief-5000' in Layer at line 24 of
'/usr/share/osm-styles/OpenTopoMap/mapnik/opentopomap.xml'

I have installed the packages gdal-bin python-gdal libgdal-dev. Is there any special configuration for gdal datasources required or a possibility to get more fine grained log-messages to estimate the reason, why this datasource file cannot be accessed?


texts: How to modify the "Network Analysis Library" in QGIS?

texts:
I want to modify the network analysis library because th graph builder does something I don't want. When there is an intersection between to lines, it creates a node in it. I want to avoid this behaviour, so I did something like :

Calculate the coordinates of every beginning and end of a line in my .shp file

Create a shp that has every node as a point

Associate the first shp containg the lines the id of the beginning point and end point in new attributes called "from" and "to"

define a new function that will each time I want to calculate the graph. It initialises a new Qgsgrapg(), then I loop over the node shp and use the addVertex() method, when it ends, I loop over my line shp and use the addEdge() method. in the end I return the graph.

My problem is it that this last function consumes a bit of time. I got over something like 300k lines, and I guess there is some python overhead in the loop. (I don't know much about these kind of things, I just assume from what I read)

So I thought that now that I made a working project, I could try to optimize it. I'm trying to understand how this function works :

http://qgis.org/api/classQgsLineVectorLayerDirector.html#af14ff34e31a2860edd156090c0f746ef

I suspect it's the one responsible of my problem. I don't know how it really works, but I guess I should find the lines responsible of the new nodes in intersections, delete them (maybe it would need some code adjustment), find a way to recompile this (There are some small c++ plugin guides in github), and tell my module from where to import the new dll (I guess that compiling this would lead to a .dll, I don't know much about c++, I have enough troubles to write std::string)

Ps: I tried to use the networkx library and it's read_shp() method to see if it's fast enough, I got a memory after 2min of processing.


texts:
I'm trying to clip a raster using a polygon but I end up with a result that doesn't seem right and where some raster pixels are outside the polygon used for the clip.

See image below where polygon is the green outline and each pixels centroid is red:

I'm using the following SQL statement:

select st_union(st_clip(a.rast,b.geom, true))
from temp_jonyp2.dhm a,
    (
    select wkb_geometry geom
    from i.gdk_bygning
    where id=347910
    ) b
where st_intersects(a.rast,b.geom)

I'm running PostgreSQL 9.6 and PostGIS 2.4.

Output of postgis_full_version is:

POSTGIS="2.4.0 r15853" PGSQL="96" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.4, released 2016/01/25" LIBXML="2.9.1" LIBJSON="0.11" RASTER

This question might relate to PostGIS polygon clip inaccurate using ST_Clip and ST_Union from tiled raster, but I can't figure out if the issue has been resolved for the version of PostGIS I'm running.

How do I clip a raster with a polygon so that only pixels inside the polygon are included?

UPDATE:
From Vince's suggestion in the comments I tried with a negative buffer. This did shave of the pixels touching the polygon in the areas where the clip looks right, but didn't help with the "wrong" parts.
See image:


texts: Mean values from a table in Earth Engine

texts:
I have uploaded a table containing information on river basins within countries (https://code.earthengine.google.com/?asset=users/basins/TFDD_BCU). I want to calculate the mean population density in the basins from this image:

var Pop1975 = ee.Image('JRC/GHSL/P2016/POP_GPW_GLOBE_V1/1975');

I have tried this code:

var BCU = ee.FeatureCollection('users/basins/TFDD_BCU');

var Pop1975mean = Pop1975.reduceRegions({
  collection: BCU,
  reducer: ee.Reducer.mean(),
  scale: 500,
});

Export.table.toDrive({
  collection: ee.FeatureCollection([
  ee.Feature(null, Pop1975mean)
  ]),
  description: 'Pop1975mean',
  fileFormat: 'CSV'
});

When running the export to the drive the process is stopped after a short while, reporting an 'internal error' but no additional information. Is the dataset too large or is there a mistake in the code?


texts: Extracting average cell values of multiple buffer zones in ArcMap?

texts:
I am trying to find a tool capable of extracting the average values of multiple buffers zones that I created around some GPS coordinates.

I have a layer in ArcMap that shows temperature in a field, I also have created a 12-ft buffer zone around some random locations (inside that field). I am aware that in order to extract multiple values in ArcMap, I can use the "Extract Multi Values to Points" (Spatial Analyst)tool; however, what I need is to extract the average values inside my multiple buffer zones.

Does anyone have any idea on how to accomplish this using ArcMap?


texts:
I am trying to run a processing script without a user defined layer output. The example below should create a vector layer with one point. The script runs without any error but the layertree is empty. If you run an other script afterwards, for example "Buffer" the layer is available and you can buffer the result from the script below.

What do I miss to make the created layer available for the layertree? My guess is that QgsProject.instance().addMapLayer(vl) does not work with the alg class? I know that the correct usage of the processing script would be using a QgsFeatureSink, but I am trying to do it the way below to get more control over the result layer.

from qgis.processing import alg
from qgis.PyQt.QtCore import QVariant
from qgis.core import QgsProject, QgsFeature, QgsField, QgsGeometry, QgsPointXY, QgsVectorLayer

@alg(name="test123", label=alg.tr("test123"), group="geophysics", group_label=alg.tr("Geophysics"))
@alg.output(type=alg.NUMBER, name="OUTPUT", label='Number of features processed')
def test123(instance, parameters, context, feedback, inputs):
    """
    tries to add a vectorlayer through a processing script
    """
    vl = QgsVectorLayer("Point", "Hello", "memory")
    pr = vl.dataProvider()
    pr.addAttributes([QgsField("name", QVariant.String),
              QgsField("field1",  QVariant.Int),
              QgsField("field2", QVariant.Double)])
    vl.updateFields()

    f = QgsFeature()
    f.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(20,20)))
    f.setAttributes(["Hello QGIS", 100, 10.1])
    pr.addFeature(f)
    vl.updateExtents()
    QgsProject.instance().addMapLayer(vl)
    if vl.isValid():
       QgsProject.instance().addMapLayer(vl)
       out = 1
    else:
       out = 0
    return {"OUTPUT":out}


texts: module 'qgis.processing' has no attribute 'run'

texts:
I am using Python for my postgraduate studies and I am having
trouble running my code when I try to use processing in Qgis using the command-line interface.
I think the error is really easy to fix it but right now I do not know how to fix it. When I run my code I have this
Traceback (most recent call last):
  File "urban_canopy_element.py", line 100, in <module>
    output=general.run('qgis:clip',{'INPUT': input_layer, \
NameError: name 'run' is not defined

I do not know where the error is. Attach the source code
import os
import sys
#from os import path
from qgis.core import *
from qgis.core import (QgsApplication,
QgsDataSourceUri,QgsCategorizedSymbolRenderer,
QgsClassificationRange, QgsPointXY,
QgsProject,QgsExpression,QgsField,QgsFields,QgsFeature,
QgsFeatureRequest,QgsFeatureRenderer,QgsGeometry,QgsGraduatedSymbolRenderer,
QgsMarkerSymbol,QgsMessageLog,QgsRectangle,QgsRendererCategory,
QgsRendererRange,QgsSymbol,QgsVectorDataProvider,QgsVectorLayer,
QgsVectorFileWriter,QgsWkbTypes,QgsSpatialIndex)
####################################33
QgsApplication.setPrefixPath("/usr", False)
#import processing

sys.path.append('~/.local/share/QGIS/QGIS3/profiles/default/python/')
#################################################
os.chdir("~/rsync/avaliar_umep/")
input_file="~/rsync/avaliar_umep/Height-rec-wg.shp"
clip_polygons_file="~rsync/avaliar_umep/grilla01.shp"
output_folder="~lixo/"
input_layer = QgsVectorLayer(input_file, "input file", "ogr")
QgsProject.instance().addMapLayer(input_layer)
tile_layer  = QgsVectorLayer(clip_polygons_file, "clip_polys", "ogr")
QgsProject.instance().addMapLayer(tile_layer)
tile_layer_dp=input_layer.dataProvider()
print(tile_layer_dp)
EPSG_code=int(tile_layer_dp.crs().authid().split(":")[1])
tile_no=0
clipping_polygons = tile_layer.getFeatures()
heightshapes=[]

for clipping_polygon in clipping_polygons:
    print("Tile no: "+str(tile_no))
    tile_no+=1
    geom = clipping_polygon.geometry()
    area=geom.area()
    print(area)

    clip_layer=QgsVectorLayer("Polygon?crs=epsg:"+str(EPSG_code)+\
    "&field=id:integer&index=yes","clip_polygon"+str(tile_no), "memory")
    clip_layer_dp = clip_layer.dataProvider()
    clip_layer.startEditing()
    clip_layer_feature = QgsFeature()
    clip_layer_feature.setGeometry(geom)
    (res, outFeats) = clip_layer_dp.addFeatures([clip_layer_feature])
    clip_layer.commitChanges()

    #print(clip_layer.extent())# linha desnecessaria neste momento

    clip_file = os.path.join(output_folder,"tile_"+str(tile_no)+".shp")
    output=processing.run('qgis:clip',{'INPUT': input_layer, \
    'OVERLAY': clip_layer, 'OUTPUT': "memory:"})
    outputlayer=output['OUTPUT']

Thanks. I am looking foward for your help.
Regards


texts:
I would like to clip some tiff files using cutline from kml file.

.\gdalwarp.exe -of GTiff -co COMPRESS=LZW -cutline boundingbox.kml source.tif result.tif

It works quite well but as a result I have files with the same size as before.

Is it possible to adjust result image size to the result content?

Here I would like to get rid of that black border.

Edit:
When I add -crop_to_cutline parameter then result is even worse:


texts: Handle string or int in server object extension

texts:
How would you also check to see if the value you are passing into into an SOE is an integer as well?

This checks for string, but returns 0 when passing number:

string searchValStr;
int searchValInt;

found = operationInput.TryGetString("searchVal", out searchValStr);
        if (!found || string.IsNullOrEmpty(searchValStr))
            // check for int
            Int32.TryParse("searchVal", out searchValInt);

How would I check for int and string using same input?


texts:
I created my plugin with Plugin Builder in QGIS. It works great.
Now I want to add more functions (http://qgis4dummies.wdfiles.com/local--files/grass-plugin/01_nuovo_mapset_grass_qgis.png like in this picture - there are many functions under GRASS plugin).
How can I do it? I tried to create another plugin, but I wasn't able to join them under one plugin.


texts: OpenLayers Plugin won't render maps correctly

texts:
In QGIS, I've created a map by adding a Tiger layer of my county, a KML of a water feature, and a geo-coded CSV of points of interest. I'm trying to add another layer using the OpenLayers Overview.

When the Overview is open in my sidebar, it correctly identifies in its small window the location that I am viewing in the main window. However, when I try to add a layer, the I instead getting a repeating Google map that just appears as a band across the middle of the screen instead of a layer.

Any ideas?


texts: Possible to add multiple fields in single arcpy statement?

texts:
I'm trying to add multiple fields to a table using arcpy and am wondering if it's possible to have multiple fields declared in the AddField management, and their values, separated by commas or something? Here is the code as stands, which is slightly cumbersome:

import arcpy, sys

arcpy.env.workspace = "C:\\gislab2\\Python\\take_home"
shapefile = "USCancer2000.shp"

#Set fieldname variables
field_name1 = "Crdrt"
field_name2 = "Crdrt_min"
field_name3 = "Crdrt_max"
field_name4 = "Ageadj"
field_name5 = "MinAgeAdj"
field_name6 = "MaxAgeAdj"

###Declare each standard population per age group, and for summed standard population
Pop1Stnd=78784000
Pop2Stnd=150486000
Pop3Stnd=45364000
TotalStnd=274634000

    #Add the fields one at a time
    arcpy.AddField_management(shapefile,field_name1,"DOUBLE",
"5","4","","","NULLABLE","NON_REQUIRED","")
    arcpy.AddField_management(shapefile,field_name2,"DOUBLE", "5","4","","","NULLABLE","NON_REQUIRED","")
    arcpy.AddField_management(shapefile,field_name3,"DOUBLE", "5","4","","","NULLABLE","NON_REQUIRED","")
    arcpy.AddField_management(shapefile,field_name4,"DOUBLE", "5","4","","","NULLABLE","NON_REQUIRED","")
    arcpy.AddField_management(shapefile,field_name5,"DOUBLE", "5","4","","","NULLABLE","NON_REQUIRED","")
    arcpy.AddField_management(shapefile,field_name6,"DOUBLE", "5","4","","","NULLABLE","NON_REQUIRED","")

I know there's a way to create a loop to do it, but was wondering if there was a simpler and more direct way.


texts: How to remove/disable a toolbar via python?

texts: How do I remove or disable a toolbar e.g., Raster, programmatically using Python?


texts: QGIS Advanced Python Field Calculator Error: File "", line 1, in  NameError: name 'DN is not defined

texts:
I am using the advanced python field calculator to change the numeric values from a field "DN" to text using the following script:
def getValue(DN):
if DN == '10':
    value = 'undifferentiated urban'
elif DN == '11':
    value = 'Urban'
elif DN == '12':
    value = 'Commercial'
elif DN == '13':
    value = 'Major Roads'
elif DN == '20':
    value = 'Undifferentiated Cropland'
elif DN == '21':
    value = 'Cropland'
elif DN == '30':
    value = 'Undifferentiated Cutblocks'
elif DN == '31':
    value = 'Graminoid Cutblock'
elif DN == '32':
    value = 'Shrubby Cutblock'
elif DN == '33':
    value = 'Tree Cutblock'
elif DN == '40':
    value = 'Undifferentiated Burn'
elif DN == '41':
    value = 'Graminoid Burn'
elif DN == '42':
    value = 'Shrubby Burn'
elif DN == '43':
    value = 'Treed Burn'
elif x == '44':
    value = 'New Burn'
elif DN == '51':
    value = 'Closed Black Spruce'
elif DN == '52':
    value = 'Closed Jack Pine'
elif DN == '53':
    value = 'Closed White Spruce'
elif DN == '54':
    value = 'Closed Undifferentiated Conifer'
elif DN == '55':
    value = 'Closed Deciduous'
elif DN == '56':
    value = 'Closed Coniferous Leading Mixedwood'
elif DN == '57':
    value = 'Closed Deciduous Leading Mixedwood'
elif DN == '58':
    value = 'Closed Mixwood'
elif DN == '61':
    value = 'Closed Shrub'
elif DN == '63':
    value = 'Closed Upland Shrub'
elif DN == '71':
    value = 'Fescue Grassland'
elif DN == '72':
    value = 'Mixed Grassland'
elif DN == '81':
    value = 'Emergent Wetlands'
elif DN == '82':
    value = 'Graminoid Wetlands'
elif DN == '83':
    value = 'Shrubby Wetlands'
elif DN == '85':
    value = 'Lichen Bog'
elif DN == '86':
    value = 'Black Spruce Bog'
elif DN == '88':
    value = 'Undifferentiated Wetland'
elif DN == '91':
    value = 'Water'
elif DN == '101':
    value = 'Permanent Ice and Snow'
elif DN == '102':
    value = 'Rock Talus'
elif DN == '103':
    value = 'Exposed Soil'
elif DN == '105':
    value = 'Upland Dune Field'
elif DN == '112':
    value = 'Cloud'
elif DN == '113':
    value = 'Haze Shadow'
elif DN == '152':
    value = 'Open Pine'
elif DN == '153':
    value = 'Open Spruce'
elif DN == '154':
    value = 'Open Undifferentiated Conifer'
elif DN == '155':
    value = 'Open Deciduous'
elif DN == '163':
    value = 'Open Shrub'
else:
    value = 'none'
return value

However, I get the following log error:
Traceback (most recent call last):
File "C:/PROGRA~1/QGIS3~1.14/apps/qgis/./python/plugins\processing\algs\qgis\FieldPyculator.py", line 176, in processAlgorithm
exec(bytecode, new_ns)
File "", line 1, in
NameError: name 'DN' is not defined
Execution failed after 0.03 seconds
Loading resulting layers
Algorithm 'Advanced Python field calculator' finished
Eventually, I would like to include this script at the tail end of a graphical modelling task I've completed.  The field name "DN" will be the same as well.


texts: Create a raster with ranges from another raster

texts: I have a raster that contains values ranging from 1-180. I am looking for a way to create a new raster where if the values of a cell in the original raster is between 1-10 the new raster cell gets the value of 1, 11-20 gets the value of 2 etc.


texts:
I'm currently trying to apply shade relief effect to dem raster layer(geotiff format, one band with float samples, min and max are around 77 and 100) using geotools. I've tried diffrent parameters setups (factor and brigthness only flag), but effects are not even noticable. I'm wondering if it's my fault or this feature isn't implemented. I don't want to use gdal because I'd like to keep dem samples unchanged. I'm attaching code sample used to generate style applied to dem layer.

    StyleFactory sf = CommonFactoryFinder.getStyleFactory(null);
    FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2(null);
    ContrastEnhancement ce = sf.contrastEnhancement(ff.literal(1.0), ContrastMethod.EXPONENTIAL);
    SelectedChannelType sct = sf.createSelectedChannelType(String.valueOf(1), ce);
    RasterSymbolizer sym = sf.getDefaultRasterSymbolizer();
    ChannelSelection sel = sf.channelSelection(sct);

    sym.setChannelSelection(sel);
    sym.setColorMap(cm.build());
    sym.setShadedRelief(new ShadedReliefBuilder().factor(22).brightnessOnly(true).build());

    Style style = SLD.wrapSymbolizers(sym);


texts: What is a "single grid stack", relating to raster data?

texts:
I created a boosted regression tree model using data points from a DEM, TPI, and land cover raster. However, I'm having problems finding resources that describe how to fit my trained model over those rasters. I found a few examples of this using R but my model was created in Python.
I found this website: http://www.wvview.org/open_source_gis/site_renders/Spatial_ML/site/index.html#predict-to-raster-data
This is actually what I'm trying to do but I don't understand what is meant by "I have generally found that it is easier to feed in the predictor variables as a single grid stack, which is prepared using GIS software, to avoid any issues with cell misalignment or different numbers of rows and columns. Python is generally not robust to dealing with these issues."
How does one "stack" rasters together to make a single grid stack which can be imported into Python?


texts:
I can add static html code to a Leaflet popup window but can not load a Dygraph chart into a popup window. How does one load dynamic content into a Leaflet popup window?
Here is a simple jsfiddle with a graph and map.
https://jsfiddle.net/scottgeowork/pybvvo27/1/


texts: Taking data from related table and populating feature attribute table in ArcGIS Desktop?

texts:
Our forester wants trees symbolized red if they haven't been trimmed, green if they've been trimmed.  Trimming info is maintained in the Tree_Maint table, which is a one to many related table to our Trees feature class in an enterprise gdb.  You cannot symbolize a feature from a related table using ArcMap.  I would like to take records from the Tree_Maint table's "Activity" field with the value "Trimmed" for the current year using the "Activity_Date" field and populate the Trimmed records into the Trees attribute table field called "TRIM_STATUS".  This would then allow me to symbolize the trees.

The end result of this work done in ArcMap 10.4.1 will be a feature service that our forester will use via Collector to select a tree and edit the related maintenance table, which I intend to use to update the "TRIM_STATUS".

I'm assuming this might require python, of which I am a novice.


texts: Automatically change symbol modified in style manager using QGIS3

texts:
I used Style Manager to add some symbols.
Then I assigned them to a feature/layer.

But if I change the Symbol / style (only via style editor), the layer will not change automatically.

How can I achieve this?


texts: Geoserver SLD for min-max scale range display

texts:
I am new to GeoServer.

I am trying to add a scale range to a layer that I am serving as a WMS. I want the layer to only be displayed between 1:30,000 and 1:1,000,000. I am using the max/min scale denominators in the SLD but am not having much luck. It works if I have just one scale denominator, but not if I introduce a second one to create a range.

Current SLD below:

      <NamedLayer>
<Name>grid</Name>
<UserStyle>
  <FeatureTypeStyle>
    <Rule>
      <MinScaleDenominator>1000000</MinScaleDenominator>
      <PolygonSymbolizer>
        <Stroke>
          <CssParameter name="stroke">#000000</CssParameter>
          <CssParameter name="stroke-width">0.2</CssParameter>
        </Stroke>
      </PolygonSymbolizer>
    </Rule>
    <Rule>
      <MaxScaleDenominator>300000</MaxScaleDenominator>
      <PolygonSymbolizer>
        <Stroke>
          <CssParameter name="stroke">#e32222</CssParameter>
          <CssParameter name="stroke-width">0.2</CssParameter>
        </Stroke>
      </PolygonSymbolizer>
    </Rule>
   </FeatureTypeStyle>
</UserStyle>

This code currently does the reverse of what I am trying to achieve. The layer is displayed beyond the two scales not between them.


texts: ArcPy(10.1) distanceTo using dynamic geometries?

texts:
I am struggling to get the arcpy point class method distanceTo to function with dynamically generated geometries with a Python addin.

I have used the default tool addin function and set self.shape (in __init__) to shape.  Selecting a polyline feature class in a drop down enables it.

class DrawLine(object):
    """Implementation for testaddin_addin.drawline (Tool)"""
    def __init__(self):
        self.enabled = False
        self.shape = "Line"

The user is able to draw a line, with a spatial reference automatically set to the data frame SR.  I use the line to perform a SelectLayerByLocation query and extract the intersecting geometries to a feature layer (MakeFeatureLayer).

For the tool, the order of intersection is important, so I grab the starting point of the drawn line and create a new point geometry:

start = arcpy.PointGeometry(arcpy.Point(line_geometry.firstPoint.X, line_geometry.firstPoint.Y))

Finally I iterate over the intersecting features and attempt to call distanceTo:

    while row:
        geom = row.getValue(shpfieldname)
        print type(start)
        print start.distanceTo(geom)
        row = rows.Next()

start is type <class 'arcpy.arcobjects.geometries.PointGeometry'>, so I know that I am working with a point geometry.

I am returning:

 File "c:\arcgis\desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py", line 809, in distanceTo
    return convertArcObjectToPythonObject(self._arc_object.DistanceTo(*gp_fixargs((other,))))
ValueError: <geoprocessing describe geometry object object at 0x2868B4A0>

Has anyone gotten distanceTo to work?  How about with a dynamic geometry?

Note: The input shapefile is simple polylines, single part, composed of pairs of vertices, i.e. as simple as possible.


texts: Solving topological errors like invalid geometries in QGIS?

texts:
I use QGIS, version 2.18.9.
I´ve several polygon-shapes with topological errors. Before I tried to apply the "intersection"-tool and got the message: Input layer A contains invalid geometries (feature 77). Unable to complete intersection algorithm.
Firstly, I used the "check-validity"-tool, which detected some line-intersection and double nodes.

Then, I tried the "topology checker" with "must not have invalid geometries".

I´ve had the impression, that all errors where found then. So I applied the LWGEOM algorithms "Make valid" but got a error another messages ([error 87] falscher Parameter see log for more details).

I´m an absolute QGIS-beginner. So, I have not a single idea how to solve this problem. It has to be a simple automatically way, ´cause I´ve a big dataset here.


texts: Combining dynamic service with map cache using ArcGIS Server

texts:
What I would like is a single ArcGIS Service that caches at the largest levels (state, regional, city) but is then just dynamic at sub city levels. From what I'm reading and understanding -- I think this is not possible -- one would have to have two separate services that display at different levels.

And yes, I know you can do the "create tiles on demand", but space, speed issues and the fact that I can't even get it to work at the moment is not ideal.

Is this possible?


texts: Improving ArcMap Batch Production of PDF Maps Performance?

texts:
I am producing 1,500 PDF maps and was wondering what changes I can make to either my machine or to ArcMap to help speed things along. Currently, it takes 5:45 - 6 minutes to produce four maps. If I change the output format to PNG, this goes down to 1:46 minutes, but I need to output PDFs.

My specs are:
-one geodatabase
-20 or so different feature classes
-python script to cycle through 30 different layers for all 50 states
-running on local machine
-ArcMap 10.6.1

What I've tried:
-clearing cache in ArcMap
-rebooting machine
-enabling hardware acceleration and background processing
-attribute indexes on query fields
-spatial indexes added
-ArcGIS cache manager declared high priority in Windows Task Manager

Any suggestions/tips that I can try to speed this up?


texts: 'st_touches' from 'sf' R package does not return all boundaries

texts:
I am trying to get the total length of all the shared boundaries between cells in a grid (following this post). However, I noticed that some of the objects have no associated boundaries, even though they do exist.

Because the issue might be related to my specific objects, you can download the initial object (study.area) here. The code bellow allows you to create the grid and the boundaries between cells.

require(sf)
require(ggplot2)

load("file_path/study_area.RData") #load study area into R, repace "file_path" by your directory

# create hexagonal grid within study_area
cell_diameter <- sqrt(2 * 10000 / sqrt(3)) # set grid area to 10000m (1ha)

grid <- st_make_grid(study.area, cellsize = cell_diameter, square = FALSE) # create grid using "st_make grid" in sf package

grid <- st_sf(grid) # convert to sf

grid$id <- 1:length(grid$grid) # add na indentifying id to each grid cell

Touching_List <- st_touches(grid) # create list identifying all the neighbouring grids

By calling Touching_List you can see that, for example, list object 1, 2 and 3 are empty, that said, they have no boundaries.

Touching_List
    Sparse geometry binary predicate list of length 16480, where the predicate was `touches'
first 10 elements:
 1: (empty)
 2: (empty)
 3: (empty)
 4: 6
 5: 7
 6: 4
 7: 5
 8: 9
 9: 8
 10: (empty)

So, I decided to have a better look on this by plotting grid cell 3 and its neighbours.

  # build demonstration plot
  ggplot(data = grid[c(3, 1, 6, 4),]) +
  geom_sf(fill = c("orangered", "peachpuff1", "peachpuff1", "peachpuff1")) +
  geom_sf(data = shp.avz, fill = NA, size = 1, colour = "black") +
  geom_sf_label(aes(label =  id)) +
  coord_sf(xlim = c(-27570, -27325), ylim = c(6935, 7260), expand = FALSE) +
  theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())

As you can see in the picture, cell 3 does have neighbours although this is not shown in Touching_List. Finally, the last step is to check whether grid cell 3 actually appears as a neighbour on cell 6, 4 and 1.

Touching_List[c(1,4,6)]
[[1]]
integer(0)

[[2]]
[1] 6

[[3]]
[1] 4

And it does not. Finally, I noticed that this only happens with cells located at the limit of my study area (see Figure 1).

EDIT 1

This Edit uses st_precision as suggested in the answer provided by Edzer Pebesma. Although setting precision to 0.1 does seem to retrieve the expected results, setting it for instance to 0.01 returns unexpected boundaries.

# Precision set to 0.1
st_precision(grid) <- 0.1
Touching_List <- st_touches(grid)
head(Touching_List)
[[1]]
[1] 2 3 4

[[2]]
[1] 1 4 5

[[3]]
[1] 1 4 6

[[4]]
[1] 1 2 3 5 6 7

[[5]]
[1] 2 4 7 8

[[6]]
[1]  3  4  7 10

As you can see in Figure 1, st_touches now returns the expected neighbours for cell 3 (cells 1, 4 and 6).

Now, let's instead set precision to 0.01.

 # Precision set to 0.01
st_precision(grid) <- 0.01
Touching_List <- st_touches(grid)
head(Touching_List)
[[1]]
[1] 2 3 4 6

[[2]]
[1] 1 4 5 7

[[3]]
[1] 1 4 6

[[4]]
[1]  1  2  3  5  6  7 10

[[5]]
[1]  2  4  7  8 11

[[6]]
[1]  1  3  4  7 10 15

In this case, we get more neighbours than expected. Let us take cell 6 as an example, as you see in the previous list st_touches returns cells 1, 3, 4, 7, 10 and 15 as neighbours. Plotting this cells

ggplot(data = grid[c(6, 1, 3, 4, 7, 10, 15),]) +
  geom_sf(fill = c("orangered", "yellow", "peachpuff1", "peachpuff1", "peachpuff1", "peachpuff1", "yellow")) +
  geom_sf(data = shp.avz, fill = NA, size = 1, colour = "black") +
  # geom_sf(data = shp.avz, aes()) +
  geom_sf_label(aes(label =  id)) +
  coord_sf(xlim = c(-27570, -27200), ylim = c(6935, 7450), expand = FALSE) +
  theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())

From the picture, cell 1 and 15 are clearly not neighbours of cell 6, although they are identified by st_touches when precision is set to 0.01. Any thoughts on why this happens?


texts: project hdf file with geospatial_lat and lon in batch process

texts: Is there any way that I can do a batch process in projecting hdf files wherein I could get and include only the geospatial_lat and lon in the batch script since every hdf file have different geospatial extent.


texts: Nonlinear regression in Google Earth Engine on each pixel in an image collection using GEE's servers

texts:
I am attempting to optimize the parameters for a double-logistic function on an annual MODIS NDVI time series in Google Earth Engine. I have previously looked at Non Linear Regression in Google Earth Engine on FeatureCollection, which is a similar problem without much clarity, other than hinting that it is possible to optimize the parameters for nonlinear regression within GEE. Previously, I have been able to use Python to change each image of the collection into an array, stack the arrays, then subset each length 46 array from each pixel and fit my function using lmfit.

However, I run into memory issues when attempting to stack arrays when using it on more than the 1x1 degree block of images I initially got it working on. So, my thought is that I need to do this without manipulating the Image Collection to a stack of np.arrays as that then attempts to load the Image Collection into local memory.

-Is this only a problem that can be solved by cutting the arrays into chunks in Python, or can I leverage GEE to reduce my computing times, beginning with the script below?

//Date range for image collection
var date_start = ee.Date('2002-01-01');
var end_date = ee.Date('2002-12-31');

//Region of interest
var roi = ee.Geometry.Rectangle(-100.1, 40, -100, 40.1);

//Clipping function for image collection
var clipping = function(image){
  return image.clip(roi)};

// define the image collection
var collection = ee.ImageCollection('MODIS/006/MOD09Q1')
    .filterDate(date_start, end_date)
    .map(clipping);

// NDVI function
var NDVIcalc = function(image){
  var red = image.select('sur_refl_b01');
  var nir = image.select('sur_refl_b02');
  var ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI');
  return ndvi;
};

//apply NDVI function to collection
var series = collection.map(NDVIcalc);

//6 parameter double logistic fundtion from Beck et al. 2006 where; t = doy, wNDVI = Annual minimum NDVI,
//wNDVI = Annual maximum NDVI, mS = Spring slope, S = SOS, mA = Fall slope, A = EOS
var dbl_logistic_model = function(t, wNDVI,mNDVI,mS,S,mA,A){
    return wNDVI + (mNDVI - wNDVI)* ( 1/(1+Math.exp(-mS*(t-S))) + 1/(1+Math.exp(mA*(t-A))) - 1)};

//Function to create range of numbers
var Arange= function(a, b, step){
    var A= [];
        A[0]= a;
        step= step || 1;
        while(a+step<= b){
            A[A.length]= a+= step;
        }
    return A;
};

//Initial parameter guesses
var params = [0.1,0.7,0.1,100,0.1,250]
//create day of year variable for "t" in model
var xdata = Arange(1,361,8);

//In python here for a length 46 array within a single cell the normal process here I would take is to wrap dbl_logistic_function in lmfit.Model and use the .fit function from lmfit to optimize parameters
//Where somefunction uses an optimizer like Levenberg-Marquardt to optimize parameters for each pixel, and mod_fit is the 6 parameter output
var mod_fit = somefunction(series, x=xdata, params)


texts:
I am working on a project to identify weeds.  I used the grid to separate the area.  Each cell is marked weed or no weed.
I would like to get the GPS of the center of each cell with a weed in it.
Does anyone have any idea how to do this?
We ultimately hope to upload the coordinates to a drone and have it spray each marked grid.


texts: Finding overlapping polygons of two Spatial Polygon Data Frames using R

texts:
I am new to GIS and GIS in R. I have two Spatial Polygon Data Frames (spdf), one contains the marine ecoregions of the world (MEOW), the other one contains Exclusive Economic Zones (EEZ) of the world. I want to find which polygons of this two spdf's overlaps (I am not sure about the term, what I mean is, if two polygons from the two spdf's share a common area, I want to find the names of those two polygons).
I imported two spdf's with readOGR function, then I tried st_intersects() and gIntersection() functions but it did not work for me.
How can I reach my goal?
Is there a function or package that handles a situation such as this?


texts: What is better for spatial data analysis? REST or SOAP?

texts:
I am working on agriculture forecasting project. It is the integration of SOA, web services, ontologies etc. What is better for spatial data analysis ? REST or SOAP ? WMS, WFS are used but either use REST or SOAP ? PLZ reply as soon as possible.

I am going to use QGIS. So please give me suggestion in that context.

Thank You.


texts:
I am trying to limit the list of available CRS for the WMS GetCapabilities of a single workspace to reduce the size of resulting XML for parsing on the front-end app side. I have found the Limited SRS list option in the WMS admin interface, added the following: 4326, 3857

This makes the change valid when going to: ...geoserver/wms?request=GetCapabilities&...

result is:

<!-- Limited list of EPSG projections: -->
<CRS>EPSG:4326</CRS>
<CRS>EPSG:3857</CRS>
<CRS>CRS:84</CRS>

But both GetCapabilities of a single workspace or a single layer do not respect that.

Going to ...geoserver/wstest/test/wms?request=GetCapabilities&...
I can really see the two CRS as the options of that particular layer,

<Layer queryable="1" opaque="1">
    <Name>test</Name>
    <Title>test</Title>
    <Abstract/>
    <KeywordList>
        <Keyword>test</Keyword>
        <Keyword>WCS</Keyword>
        <Keyword>ImageMosaic</Keyword>
    </KeywordList>
    <CRS>EPSG:3857</CRS>
    <CRS>CRS:84</CRS>

but I still also receive the basic full list of the WMS service itself.

<Layer>
    <Title>GeoServer Web Map Service</Title>
    <Abstract> A compliant  implementation of WMS plus most of the SLD extension (dynamic styling). Can also generate PDF, SVG, KML, GeoRSS </Abstract>
    <!-- All supported EPSG projections: -->
    <CRS>AUTO:42001</CRS>
    <CRS>AUTO:42002</CRS>
    <CRS>AUTO:42003</CRS>
    <CRS>AUTO:42004</CRS>

Etc...

Should this change not be sufficient to remove the additional CRS from the resulting XML, or is there a way to propagate this change to the workspace/layer GetCapabilities response separately that I missed?


texts: Google Earth Pro does not completely open TAB files created in MapInfo 15+

texts:
From time-to-time, I have the need to open tab files I have created in MapInfo with Google Earth Pro.

Aside from minor hiccups, I have never had an issue with importing a Tab file.  The process I always follow in Google Earth Pro is:

  File -> Open -> Select *.tab from File Type -> Select the file and pick the style template.

Recently, however, I have upgraded to MapInfo 15 (MapInfo 15.2 64-bit to be precise), and since then whenever I try to open a tab file in Google Earth, I only get the first item or 2 from the tab file and the remaining features (usually polygons) are lost.  If I look at the tab in MapInfo, the data is present, so it seems MapInfo 15 is doing something differently that causes the tab file to not completely import into Google Earth Pro (GE version 7.1.5.1557).

Reverting to MapInfo 12 is not really an option, so is there anything else I can do to allow for Google Earth to completely parse these files?

To provide some additional details.  I can take the MapInfo 15.2 created tab file and open it in MapInfo 12 without an issue on another machine.  If I do "Save Copy As" from MapInfo 12 and create a new file, I can then open that file in Google Earth Pro just fine.  But anytime I directly open a MapInfo 15.2 created tab file in Google Earth, I only get the first entry, sometimes 2 or 3.  The remaining data is lost.

To add some additional info that I have come across recently, without upgrading either piece of software (still at the versions identified initially), I tried to open a list of locations identified by Lat/Lon that was created by MapInfo 15.  There are 153 locations in the TAB file, but when I imported it into Google Earth, I got an error message from Google Earth saying:

  Some features were not important due to invalid or out of range coordinates

and I only got 35 locations of the 153.  All of the coordinates appear to be valid and there is only maximum span of just over 3 degrees for both Latitude and Longitude, so it seems illogical that there is actually something invalid in the TAB file (opening in MapInfo 12 and resaving it allows it to be imported to GE correctly without issue).

This is an example of a MapInfo 15.2 tab file that does not import completely into Google Earth

!table
!version 950
!charset WindowsLatin1

Definition Table
  Type NATIVE Charset "WindowsLatin1"
  Fields 2
    Column1 Float ;
    Column2 Float Index 1 ;
begin_metadata
"\MapInfo" = ""
"\MapInfo\TableID" = "31fb8af6-d29c-4c9c-a8ec-bf9e82e3b343"
"\MapInfo\ParentTableID" = "644f95b8-3201-4436-bdce-bdcff25c2a11"
end_metadata

If I take a file that I have opened and resaved via MapInfo 12, the contents of the TAB file are identical, even though it will correctly open in Google Earth

I'm not an expert in the structure of the MapInfo file formats, but this leads me to believe that the issue is in the DAT, ID, or MAP files. It's almost like Google hits the end of a line and thinks it is the end of the file.  Was there any structure changes in the MAP or DAT files that could explain this behavior?


texts: Creating CRS of ISG (superseded NSW coordinate system) in QGIS?

texts:
I want to use a superseded projection (NSW integrated survey datum - referred to as ISG) because my main client utilises it. The projection is based on AGD66. We regularly receive survey files in this projection and need to prepare designs in this projection.

Unfortunately no one I know, including surveyors,  seems to know much about the ISG projection, let alone set up a CRS for it.

When I import ISG data into an existing QGIS file, the files reproject fine and appear in the right location.  However the coordinates are obviously incorrect, making it difficult to work with.

Is there a way to extract and set a new datum based on data I know is already in ISG format?


texts:
Fig.1 is the image in zoom of 15.5% ,and Fig.2 is fixed zoom in 59.0%.

Obviously,Fig.1 window fails to display.However,the image can display regularly after zooming into a certain percentage.

What is the reason for display failure of Fig.1?

                                        Fig.1

                                        Fig.2


texts: r - sf::st_cast("LINESTRING") - keeping first linestring only warning

texts:
I have a dataset with geometry column in which there are linestrings and multilinestrings. While keeping the linestrings I want to convert the multilinestrings to linestrings (which should potentially increase the number of rows of the sf dataframe). Unfortunately when I use sf::st_cast("LINESTRING") there is a warning telling me that it is getting rid of all except the first linestring when transforming. Is there a way to keep all linestrings from the multilinestring when using sf::st_cast. Reproducible example with warning below:

library(sf)
library(dplyr)

# sample dataframe - creating linestrings
df1 <- data.frame(lon = 1:10, lat = 1:10, var = c(1,1,1,2,2,2,3,3,4,4)) %>%
  st_as_sf(coords = c("lon", "lat"), dim = "XY") %>% group_by(var) %>%
  summarise(geometry = st_union(geometry), do_union = F) %>%
st_cast("LINESTRING")

# creating a multilinestring
df2 <- df1[1:2,] %>% mutate(var = c(1,1)) %>% group_by(var) %>%
  summarise(geometry = st_union(geometry), do_union = F) %>%
  st_cast("MULTILINESTRING")

# combining the two
df <- rbind(df1, df2)

# trying to convert only the multilinestring to two linestrings not changing
the already existing linestrings
df <- df %>% st_cast("LINESTRING")

# Warning message:
# In st_cast.MULTILINESTRING(X[[i]], ...) : keeping first linestring only

I can do it manually first converting everything to multilinestring, and after everything to linestring like in the following:

df <- df %>% st_cast("MULTILINESTRING") %>% st_cast("LINESTRING")

but is there maybe a better way of doing this?


texts: Is there a way to convert multiple columns at once from string to text, or to other data formats?

texts:
I have multiple data fields that need to be converted from string to text files. When I save them as "text" in excel or open office, they still load up as string files in QGIS. How can I convert more than one column at once from string to text (or even integer, etc. if the situation arises)?

Thanks for your help!
Alex


texts: Multiple histogram from a raster image in R

texts:
I need to get multiple histograms from a  tif image. I have the following data:

A big raster tif image:

tif1 <- raster("myRaster.tif")
hist(tif1, breaks= 10) # gives the entire histogram

A shapefile of a country. This shapefile contains a polygon for each administrative area:

shp_np <- readOGR("NPL_adm2.shp")

I can think of a way to get the histogram by repetitively cropping raster tif using each Polygon in shapefile one by one. But, it is not efficient.

Now, I am searching for easy way to get the histogram information for each polygon in the shapefile.


texts: Does Leaflet JS API provide a Geocoding service?

texts: I am considering replacement of Google Maps with Open Street Maps and the Leaflet JS API. Does the Leaflet JS API provide a Geocoding service?


texts: Stack sat images and run operation on all

texts:
I have about 300 MODIS satellite with different acqusition dates, covering the same area.
I am not used to work with so many tiles at once so I ask for advice. I am working with ArcMap10.1

What I want to do is:

Stack all the images on top of each other
clip them with the same extent.
Run raster calculator for each of them, with the same function.

I guess my biggest problem is how to stack them in a good way so I later can run my tools on the whole pack.

Best regards!


texts:
I tried compositing by two weeks.
var modis = ee.ImageCollection('MODIS/006/MCD43A4')
  .select('Nadir_Reflectance_Band2');

var twoweeks = ee.List.sequence(1,365,14);
print (twoweeks);

var byTwoWeeks = ee.ImageCollection.fromImages(
      twoweeks.map(function(w) {
        return modis.filter(ee.Filter.calendarRange(w, w, 'DAY_OF_YEAR'))
                    .mean().multiply(0.0001)
                    .set('week', w);
}));

print (byTwoWeeks);

Output would be 27 items, one for every two-weeks.
Is this code correct?
How can I check that images are grouped correctly?
When I tried as follows.
var byTwoWeeks = ee.ImageCollection.fromImages(
      twoweeks.map(function(doy) {
        return modis.filter(ee.Filter.calendarRange(doy, doy, 'DAY_OF_YEAR'));
}));

print (byTwoWeeks);

ImageCollection (Error)
ImageCollection.fromImages, argument 'images': Invalid type. Expected type: List<Image<unknown bands>>. Actual type: List<ImageCollection>. Actual value: [<ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>]


texts: Getting all web services and database connections from QGIS in PyQGIS?

texts:
I need to write a plugin for QGIS but I am not very familiar with QGIS itself and kind of new to all of this.

I want to get all connections programmatically. By connections I mean everything from WMS to WFS to PostGIS.

I have discovered that those connections are stored in the QGIS3.ini but have not yet found a way in Python to get this data (except for scanning trough the file for the information, but surely there most be a better way).


texts: Cropping or removing borders (not NoData) from multiple scanned and georeferenced maps (ArcMap or QGIS)

texts:
I have several hundred scanned maps in TIFF format, each correctly georeferenced and forming a large grid of images. The problem is that the border for each map is present - a white band, with a black edge (not NoData), containing references etc, and overlaps the adjacent tiles. I want to remove these borders from all maps (see below - the white cross through the middle and the edges around the outside).

Is there a way to remove X number of pixels from the top, left, bottom, and right of rasters? And is there a way to do it for multiple rasters in the same process? I.e., rather than using the extent (as with the Clipper tool). It would be very easy in an image processing program like IrfranView, but I want to retain the geo-data (World files etc) which would be thrown off as the extent and position of the corners would change.

I can use QGIS or ArcMap.


texts:
I create api function for searing parcels and create call from SearchProviders. All it's working fine when i try to search, but when I select object form searchbox its only zoom map to the parcel but don't highlight the parcel polygon. Also in qgis project I checked the option Add geometry to feature response
This is my searchprovider
function ParcelLocationSearch(text, requestId, searchOptions, dispatch) {
    axios.get("http://127.0.0.1:3003/v1/query/parcel?columns=id%2Cnum_parcel%2Ckat_opst_1%2CST_Extent(geom)%20as%20bbox&filter=num_parcel%20like%20'%25"+text+"%25'&sort=num_parcel&group=id%2Cnum_parcel%2Ckat_opst_1" )
    .then(response => dispatch(ParcelLocationSearchResults(response.data, requestId)));
}

function ParcelResultGeometry(resultItem, callback)
{

     axios.get("http://127.0.0.1:3003/v1/query/parcel?columns=ST_AsText(geom)%20as%20geom&filter=id%3D"+text)
    .then(response => callback(resultItem, response.data, "EPSG:4326"));
}

function ParcelLocationSearchResults(obj, requestId)
{
    let results = [];
    let currentgroup = null;
    let resultGroups = {};
    let groupcounter = 0;
    let counter = 0;

    (obj || []).map(entry => {

        resultString = (entry.bbox.replace('BOX(','').replace(')','').replace(' ',',').replace(' ',',')).split(',');

        let x = 0.5 * (parseFloat(resultString[0]) + parseFloat(resultString[2]));
        let y = 0.5 * (parseFloat(resultString[1]) + parseFloat(resultString[3]));
        if(groupcounter ==0) {
        currentgroup = {
            id: "ustergroup" + (groupcounter++),
            title: entry.displaytext,
            items: []
        };
        results.push(currentgroup);
    }
      else  {

            currentgroup.items.push({
            id: entry.id,

            layer:"parcel",
            text: entry.num_parcel,
            label:entry.num_parcel,
            title:entry.num_parcel,
            marker:true,
            x: 0.5 * (parseFloat(resultString[0]) + parseFloat(resultString[2])),
            y: 0.5 * 0.5 * (parseFloat(resultString[1]) + parseFloat(resultString[3])),
            crs: "EPSG:4326",
            bbox: [x, y, x, y],
            provider: "Parcel"
        });
    }
    });

    return addSearchResults({data: results, provider: "Parcel", reqId: requestId}, true);
}
SearchProviders: {
        "coordinates": {
            labelmsgid: "search.coordinates",
            onSearch: coordinatesSearch
        },
        "Parcel": {
            label: "parcels",
            onSearch: ParcelLocationSearch,
            getResultGeometry: ParcelResultGeometry
}
}


texts:
I can't wrap my head around how to change this legend name.
The left picture shows the raster in the Toc with that annoying "VALUE", the centre shows how the legend appears and the right one shows what I see by going to layer properties.

I would just like to omit river_density and VALUE and just display a name I choose for that legend.  The raster was created with the line density tool.


texts: How to handle layers with different spatial references in ArcGIS Runtime SDK WPF?

texts:
I am currently writing the "Zoom to Layer" function for my application. However, the spatial references between my imported layer and my base map are different.

Is there a way around this? Can I somehow project all of my shapes/graphics so that they are in the base map's spatial reference?

Is my only option to remove all layers, change spatial reference, and add layers back in?

It seems like ArcEngine's spatial reference had less restriction and was easier to use. Any insights?

Thanks for all your help,
Kev


texts: What are some web based clients for GeoServer?

texts: There is a list of Web Clients for Geoserver on the Geoserver site (http://geoserver.org/display/GEOSDOC/Clients). Unfortunately, it is posted as no longer maintained.  Does anyone have a current list of web clients?


texts: Randomly move (and rotate) an sf polygon within a broader extent using R

texts:
I have an area of interest within which I have multiple polygons (two in the example below). I would like to randomly move each of the polygons so that they are completely within the larger extent and do not overlap the original placement. I am working within sf.

Using the nc simple feature, we can run the following:

library(sf)
library(dplyr)
nc <- st_read(system.file("shape/nc.shp", package="sf"))

#Make outline extent
Extent <- nc %>%
  mutate(State = "NC") %>%
  group_by(State) %>%
  summarize()

#Make two polygons
Poly <- nc %>%
  filter(NAME %in% c("Wilkes", "Greene"))

As reproducible example, I would like to randomly move both the "Wilkes" and "Greene" polygons within the Extent. It is fine if the new location of "Wilkes" overlaps the old or new (random) location of "Green" but the random location of "Wikles" should not overlap the old location of "Wikels".

The end result would be an sf multipolygon with 4 polygons, the two originals and the two random placements (ideally within the same "Name" grouping, but with a separate field delineated which is original and which is random).


texts:
I'm trying to replicate something that was fairly easy in OpenLayers 2 in OpenLayers 3, which is: dragging an existing feature to a new location on the map.

In OpenLayers 2 I just added a DragFeature control to my map like this:

this.mapRasterDragFeature = new OpenLayers.Control.DragFeature(_self.rasterVectorSource,
     {'onComplete': function(e) { console.log(e);}});
this.mapRaster.addControl(_self.mapRasterDragFeature);
this.mapRasterDragFeature.activate();

But this control seems to have disappeared in OpenLayers 3. I've searched far and wide in the documentation for the new way to this, but I can't seem to find it. I hope I'm overlooking something.

I'm adding the feature through a click event on the map like this:

this.mapRaster.on('singleclick',function(e){
    var iconFeature = new ol.Feature({
        geometry: new ol.geom.Point(e.coordinate),
    });
    iconFeature.setStyle(_self.iconStyle);
    _self.rasterVectorSource.addFeatures([iconFeature]);
});

All I've found looking like a drag interaction is ol.interaction.DragAndDrop(), but I think it does something different than just dragging and dropping a feature...

Can someone help me out on this one?


texts: QGIS 3: Customize Data Source Manager Toolbar

texts:
I have noticed that QGis version 3.2.1 does not have the "Add vector layer" button

on the Data Source Manager toolbar:

so each time I want to add a vector layer I need to go to Layer --> Add layer --> Add vector layer, which opens the Data Manager; or use the shortcut Ctrl+Shift+V.

Is there a way to customize the Data Source manager toolbar and add the "Add vector layer" button?


texts: Using ArcPy to run Create Random Raster tool executes successfully every other time, but fails with ERROR 999999 on alternate runs?

texts:
I am trying to create a python script to run several ArcMap tools in succession on a DEM. The script will ideally create a random raster (based on a normal distribution, with a mean of 0 and a standard dev. of 2), then add this to a DEM to form a DEM + random raster. The viewshed tool will then be carried out on this DEM+. I want this to loop 20 times.

This is the code I currently have:

# Program to calculate probable viewshed from a DEM using ArcMap tools
# Author: David Ely, 02/06/2016
# Status: Work In Progress
# This program must be placed in the same directory folder as the DEM

# Import the required libraries
import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension = 'Spatial'

# Set the environment to allow overwrites to data
arcpy.env.overwriteOutput = True

# Load in the DEM and store it's extent
print 'DEM file name: '
DEM = Raster(raw_input()) # User inputs the name of the DEM file
DEMextent = DEM.extent

# Enter the RMSE of the DEM and the output path for the analyses
print 'Enter RMSE of DEM: '
RMSE = raw_input()
print 'Enter output path for analyses: '
outpath = raw_input()

# define i for incremental loop and loop for analysis
i = 0
while i<20:
    print 'procedure #' + str(i)
    distribution = 'NORMAL 0 ' + RMSE
    outname = 'rand' + str(i)
    arcpy.CreateRandomRaster_management(outpath, outname, distribution, DEMextent, 5)
    randrast = Raster(outpath + '\\' + 'rand' + str(i))
    addrast = randrast + DEM
    i = i+1

before I added the code to carry out the map algebra in the while loop the create random raster tool would work every time. However, once the next two lines of code were added the createrandomraster tool fails to execute every other time (error 999999).

Does anyone know why this could be?


texts: Export DWG and TIFF with ArcMAP open in different locations in AutoCAD

texts:
I am trying to export some data from an ArcMAP I have, namely an aerial image (as a TIFF) and boundaries (as DWG).  Both files are on the same coordinate system.

I am trying to create some drawings in AutoCAD MAP.  I create a blank file and set the coordinate system to match the source data, and then import the two files.  However both data opens in drastically different locations, both of which dont match the ArcMAP?

I export both layers by right clicking->data->export and following the dialogs.

Strangely the DWG appears to be in millimetres when I attach it?

Is there a method to extract multiple layers in a batch process and force the coordinate system and units to match?


texts:
I have a road network file (openstreetmap exported file) that has various types of edges (freeway, highway, primary, residential, service, ...). Some of these edges have multiple lanes (parallel), so that each one of these lanes considered a separate edge in road network, for example:

How can I merge these lanes in one lane and construct the road network algorithmic (I know in some applications, we can select ways and merge them with right click, but I can't do this for 100000+ edges in large cities)?


texts: How to run Refactor Fields algorithm in Python console?

texts:
I noticed that the Processing algorithm, Refactor Fields has a fields mapping parameter. I don't think I've come across this and definitely not used it before from the Python Console.

It provides an example on how to use it:

processing.runalg("qgis:refactorfields",
                  input,
                  [{'name': output_field_name,
                    'type': output_field_type,
                    'length': output_field_length,
                    'precision': output_field_precision,
                    'expression': expression_based_on_input_layer
                  }],
                  output)

I tried to modify this in the console but I receive no output. Here is the typical code I used:

import processing
input = qgis.utils.iface.activeLayer()
output = "C:/Users/Me/Desktop//test.shp"

processing.runalg("qgis:refactorfields",
              input,
              [{'name': "Value",
                'type': 0,
                'length': 10,
                'precision': 10,
                'expression': "Value"
              }],
              output)

Anyone know where I am going wrong? I'm using QGIS 2.12.0-Lyon with Processing version 2.10.2.


texts: Identify the min and max x-y coordinates of an ee.Geometry (Earth Engine)

texts:
I am using the EE API in a Python Jupyter Notebook.
I am building a map based on the LSIB_SIMPLE country boundaries.
My objective is to display my information in a map and adapt the zoom to fit the country size (Singapore is way smaller than Congo for example)
I would like to compute both the longest north-south and east-west lines inside the country to adapt my zoom accordingly on this object:
country_code = CG
country = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017').filter(ee.Filter.eq('country_co', country_code))

geom = country.geometry()

Is there a way to find the 4 needed points in the geometry object (northest, southest, eastest, westest) so that I can compute distances between them ?


texts: Can you use Landsat 7 to make studies of land surface temperature, even though it has the scan line problem?

texts:
Is it recommended to use Landsat 7 to study the temperature of the surface?
If you fill the gap of the scan line errors, would that gap filled give you different results from the original image?


texts: How georeference a DXF file in QGIS 2.14?

texts:
In Illustrator I exported a PDF as DXF and inserted it in QGIS as Vector and obtained two layers as points and lines but it is misaligned and in a unknown location.

I know georeference a raster file, but not a PDF or DXF.

The PDF file is the first one (named "Caracterização Geral") here in this link: Caracterização Geral PDF file. It say "PT-TM06/ETRS89".

This is the DXF file in my Google Drive.

I can locate one point, the City Hall - 41.1298608, -8.6066119 [41°07'47.5"N 8°36'23.8"W]. Is the layer called "Câmara_Municipal" in PDF file, a cyan color point.

My project is EPSG: 3857 / WGS 84 / Pseudo Mercator

First layer is a Bing Aerial (Web > OpenLayersPlugin > BingMaps >
Bing Aerial) - CRS: EPSG 3857 / WGS 84 / Pseudo Mercator
Another layers with shapefiles lines and points captured from OSM. -    CRS: ESPG 4326 / WGS84, which is OK on top the base layer.


texts: QGIS Orfeo toolbox error when computing image statistics: "Could not find application "ComputeImageStatistics"

texts:
I am new to QGIS, and am trying to compute image statistics with the Orfeo toolbox (to use in a future classification project). I have downloaded the Orfeo toolbox and have it activated in my toolbar, as the image shows:

I am trying to compute the image statistics on an orthomosaic with three bands (red-green-blue). However, when I try to run the tool and save the output as a .xml file, I get the following errors:

There seems to be an issue with the file path. How could I correct this error and run the tool?


texts: QGIS version migration of Styled Layer Files

texts:
I have a stand alone application which reads SLD files, to draw something on a map canvas (development started under QGIS version 2.0.1). Unfortunally I have to deal with these old SLD files version QGIS 2.0.1 within a new context QGIS version 2.4.0 (security update Debian).

Is there a tool to migrate SLD files from 2.0.1 to 2.4.0?

Error or Warning?

Loading a style file that was saved with an older version of qgis (saved in 2.0.1-Dufour, loaded in 2.4.0-Chugiak). Problems may occur.

uname -a: Linux paprika 3.2.0-4-amd64 #1 SMP Debian 3.2.68-1+deb7u1 x86_64 GNU/Linux

/etc/debian_version 7.8

Sample BOR.sld:

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor xmlns="http://www.opengis.net/sld"
  xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  version="1.1.0"
  xmlns:xlink="http://www.w3.org/1999/xlink"
  xsi:schemaLocation="http://www.opengis.net/sld
  http://schemas.opengis.net/sld/1.1.0/StyledLayerDescriptor.xsd"
  xmlns:se="http://www.opengis.net/se">
  <NamedLayer>
    <se:Name>BORDER</se:Name>
    <UserStyle>
      <se:Name>BORDER</se:Name>
      <se:FeatureTypeStyle>
        <se:Rule>
          <se:Name>Single symbol</se:Name>
          <se:PolygonSymbolizer>
            <se:Stroke>
              <se:SvgParameter name="stroke">#ffff00</se:SvgParameter>
              <se:SvgParameter name="stroke-width">0.35</se:SvgParameter>
            </se:Stroke>
          </se:PolygonSymbolizer>
        </se:Rule>
      </se:FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

Sample BOR.qml

<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
<qgis version="2.0.1-Dufour"
 minimumScale="0"
 maximumScale="1e+08"
 minLabelScale="0"
 maxLabelScale="1e+08"
 hasScaleBasedVisibilityFlag="0"
 scaleBasedLabelVisibilityFlag="0">
  <renderer-v2 symbollevels="0" type="singleSymbol">
    <symbols>
      <symbol alpha="0" type="fill" name="0">
        <layer pass="0" class="SimpleFill" locked="0">
          <prop k="border_width_unit" v="MM"/>
          <prop k="color" v="0,0,0,255"/>
          <prop k="color_border" v="0,0,0,255"/>
          <prop k="offset" v="0,0"/>
          <prop k="offset_unit" v="MM"/>
          <prop k="style" v="solid"/>
          <prop k="style_border" v="solid"/>
          <prop k="width_border" v="0.8"/>
        </layer>
      </symbol>
    </symbols>
    <rotation field=""/>
    <sizescale field="" scalemethod="area"/>
  </renderer-v2>
  <customproperties>
    <property key="labeling" value="pal"/>
    <property key="labeling/addDirectionSymbol" value="false"/>
....

Is it wise to change persistent production features, like settings of a map composition, within a major release?


texts: Extracting pixels with value in one raster but NoData in other raster?

texts: I have two raster A and B, and they are representing the info for Florida. Ideally A and B should have same coverage. However, due to some data processing issue, B's boundary is a bit smaller than A's. This means there are pixels with value in A, but NoData in B. Not just along the boundary, within B's coverage, there are a few holes, meaning pixel value is NoData. My question is how I should extract these pixels to form a new raster, as well as making their pixel value same as A's.


texts: Difference between reduce(ee.Reducer.mean()) and mean() for ImageCollection in GEE

texts:
I am working with S5P data for NO2. I have an image collection and want to obtain the mean NO2 for specified filters. The following code is from the documentation.
var collection = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_NO2')

.select('NO2_column_number_density')
  .filterDate('2019-06-01', '2019-06-06');

var band_viz = {
  min: 0,
  max: 0.0002,
  palette: ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']
};

Map.addLayer(collection.mean(), band_viz, 'S5P N02');
Map.addLayer(meanImage,band_viz, 'another method');

I was wondering what is the use of reduce method on ImageCollection when there is a mean method already available.
var meanImage = collection.reduce(ee.Reducer.mean());

What is the difference between applying mean directly versus using reduce method


texts:
I have a big GeoTIFF(25GB) and I want to create tiles from it.

I know GeoServer have the feature, but I can't use it. I found GeoTools example for cropping tif files here, but it just create some tiffs from a tiff file, not creating png/jpeg tiles. I found JAI-EXT, but there is no example.
I found gdal2tiles.py and use it successfully in java, but we need to manage creating tiles(for example some times, we didn't need to create all tiles, just we want to create 20% of them in special area)

Can anyone give an example to create a special tile from GeoTIFF file?


texts: postgis, mapserver & openlayers: if data is not in table give an error / pop-up

texts:
I have created a website with postgis, mapserver & openlayers.
On the website you can enter a postal code and then a drivetime isochroon is created.

The problem now is that you can enter anything in the input box. If you enter a postal code that does not exist (or maybe some letters) nothing happens. No isochron is created, which is correct because it is not a valid code.

What I want is that if you enter a not existing postal code I would like to have a pop-up (or something like that) which says: 'this is not an existing postal code, please enter another postal code'

Someone any idea how I can create this?

The website is: maps.geodan.nl/reistijden.htm

Thanks, Stefan


texts:
I'm using ArcGIS 10.1 and developing an application using visual studio 2010 and ARCObject SDK 10.1.

I would like to use identity tool in my application and if not possible I want to execute a query on the attribute table of the map layer when a user clicks on it.

Any ideas on how I can accomplish this task? From my research,I think IQuery interface may be used; however, I would like to do the query based on a mouse click event.


texts: How is a minimal tile server implemented?

texts:
I am trying to understand how a typical tile server is implemented in practice. I am not considering setting up a tile server but rather trying to understand from a conceptual point of view the algorithms behind a typical tile server.

So for example,  as a starting point, a tile server receives a request similar to the following URL template  http://[abc].tile.openstreetmap.org/zoom/x/y.png.

What are exactly the next chain of events in the processing of such a request? In the typical case of OpenStreetMap, I know that the actual data is stored in PostGis and that mapnick is used for rendering the tile.

However I have a hard time understanding how the tile is actually rendered given the zoom, x, and y inputs.

Edit: I found the following diagram, which gives you a high-level architecture of the OpenStreetMaps tile server. http://wiki.openstreetmap.org/wiki/Component_overview


texts: Blurry menu text on Macbook Pro 15'' retina display

texts:
I'm using qgis on my macbook pro and all the text of the submenus/popup windows is displayed a little blurry/pixelated. See this example screenshot:

Is this a normal behavior on retina displays or is there a fix for this problem?

Hard facts:

macbook pro OS X El Capitan
qgis 2.12.1-Lyon (installed from http://www.kyngchaos.com/software/qgis)


texts: Customizing Create Features window in ArcMap?

texts:
Can I customize Create Features window in ArcMap? (whether by using ArcMap capabilities or by writing ArcMap-Addin)

I've read about Creating custom ArcMap TOC views, is there something similar for Create Features window?

My requirement is to visualize template names (not layer names), in addition to classify templates in tree-like organization.


texts: Python based GUI with some GIS capabilities

texts:
I would like to construct a Python based GUI which has basic GIS visualization capability. This capability should include:

1) importing OpenSteetMap

2) importing shapefiles/rasters

3) turning layers on and off

4) zooming, panning etc

I was looking into packages like qt, tkinter or mpl_toolkits. I was able to plot a shapefile using mpl_toolkits.Basemap but I have the impression it is not the most suitable tool for this kind of application... Any ideas?


texts:
In earlier versions i was able to open pdf files by using an action command like this one:

cmd /c [%path%]

this function no longer works on version 3.8. Any ideas why and what an alternative could be?


texts:
I make web application in Eclipse using GDAL. In Eclipse in Debug Configurations i add variable path=path/to/gdal/folder. In project i add gdal.jar to build path and copy gdal folder in lib directory. And all works fine. After i export my application to WAR file and want to deploy it in another computer. I download gdal for Win 64bit, put it in tomcat/shared/lib and add it in Win varible PATH.
And get error:

java.lang.UnsatisfiedLinkError: org.gdal.ogr.ogrJNI.GetDriverCount()I
org.gdal.ogr.ogrJNI.GetDriverCount(Native Method)

What i can do wrong?

UPDATE

In catalina.bat i set set JAVA_OPTS="-Djava.library.path=/tomcat/shared/lib" but its not help me. There i get answer that if i see something like this in errors its mean that library is loaded and problem in library itself.


texts: Upgrading guidelines for ArcGIS Desktop?

texts:
We have a copy of ArcGIS 10 that we are preparing to use to upgrade our existing ArcGIS 9.3 installation.

Should I uninstall ArcGIS 9.3 first?

Or do I install ArcGIS 10 over it?

Are there any best practises or guidelines I should be aware of?


texts: Labels are not appearing at certain scales

texts:
1) I wanted to let the labels of my 'river' appear in scale 25 Mio.

2) the name of the river shall not get smaller when I zoom-in and zoom-out, the size of the label should stay the same.

I want the labels to appear for 4 rivers, which are in my attribute table named as 'hoe_name'.

Now the name of the river in scale 25 Mio only one is displayed:

In scale 12 500 000 appears only three:

.

These are my label settings: 'Show labels for this layer':

How can I get all 4 (labels of my river) to appear at a scale 25 Mio. and with the same size?


texts: Openlayers interpolating a raster source

texts:
I need to show a raster layer in OpenLayers keeping its pixels as they are, i.e without iterpolating them when zooming in.

I managed to do that with an ImageStatic source following this answer: https://stackoverflow.com/a/35877192/6185152

But when using a Raster Source, the image gets interpolated. Take for example the following code that lets you show the ImageStatic layer or a Raster layer by commenting/uncommenting them (available at http://jsfiddle.net/cnadia/vfajwdmr/):

var extent = ol.proj.transformExtent(
    [-5.6342, 50.3331, 1.6607, 53.0559], 'EPSG:4326', 'EPSG:3857');

source_img = new ol.source.ImageStatic({
  url: 'https://i.imgur.com/1mDzkTJ.png',
  imageExtent: extent,
  crossOrigin: 'anonymous'
});

layer_img = new ol.layer.Image({
        source: source_img});

source_raster = new ol.source.Raster({
    sources: [layer_img],
  operationType: 'image',
  operation: function(pixels, data) {
    debugger;
    var pixel = pixels[0];
    return pixel;
  }
});

layer_raster = new ol.layer.Image({
        source: source_raster
    });

var map = new ol.Map({
  layers: [
        //layer_raster
        layer_img
  ],
  target: 'map',
  view: new ol.View({
    center: ol.extent.getCenter(extent),
    zoom: 6
  })
});

map.on('precompose', function(evt) {
  evt.context.imageSmoothingEnabled = false;
  evt.context.webkitImageSmoothingEnabled = false;
  evt.context.mozImageSmoothingEnabled = false;
  evt.context.msImageSmoothingEnabled = false;
});

what I get is:

Correct image using ImageStatic:

Incorrect (interpolated) Raster Source Layer:

I tried using the Raster Source both in 'pixel' or 'image' mode, and the result is the same.

In this example I'm "bypassing" the image (return value = input value), but I obviulsy need to modify the pixels and that is why I need the Raster Source.

Any ideas how to make this work?


texts: Where can I get data for traffic signals, stop signs, and other sign assets/inventory?

texts:
I'm working on a project involving roads.  The quality of the service would benefit from knowing if a given road is continuous or if it's broken up by stop signs, traffic signals, and other things that effect a motorist.

I haven't found anything googling around online.  I also contacted my local city government, but they only grant access to the data for city projects.  I'm guessing most areas would have similar rules.

Is there another source, public or commercial, that could help me out?  Or am I out of luck?

EDIT - I should note that I'm looking for data in the USA.


texts: Is there any tool in QGIS that traces a line like the trace tool (Editor Toolbar) in ArcGIS? I want to draw a line and trace this line back so that I get the beginning and the end of the line at the same position. Any suggestions?


texts:
I have a series of locations stored within a Vector layer and want to write a query to understand which of these locations fall within the polygons from another layer using the QGIS Python console.

So far I have got this code - however it returns all locations, some of them several times.

canvas = qgis.utils.iface.mapCanvas()
layer1 = canvas.layer(0) #layer 1 is the polygons
layer2 = canvas.layer(1) #layer 2 is the locations
featuresloc = layer2.getFeatures()
for feature in featuresloc:
    geomloc = feature.geometry()
    featurespol = layer2.getFeatures()
    for featurep in featurespol:
        if featurep.geometry().within(geomloc):
            attributes = feature.attributes()
            print attributes

Can anyone provide any pointers as to how I can correct this?


texts: Calculating area from DEM value range in QGIS

texts:
I was wondering if it is possible to calculate the area of a DEM in specific value ranges, e.g. the total area between 200 and 300 meters of elevation. Using QGIS and my DEM is 1m. Would like to do this for my DEM, so that in a second step I can visualise the number of features in each value range.
Any help is appreciated.


texts: Collapse dualline to centreline only creates outline of dataset using ArcMap

texts:
I'm trying to collapse all divided roads in one centreline and to integrate any other spurs or errors in the dataset to create one set of centrelines. I have used Merge Divided Roads but it leaves a lot to be desired (most of the roads I wanted collapsed did not merge).
Here is the dataset:

And here is what happens when it is collapsed:

I have also repaired geometry, reprojected and made it into a singlepart feature but nothing seems to work. I have tried it on multiple computers.


texts: Using geojson-vt to convert MongoDB output server side

texts:
I am using Mapbox with React.js on client-side, and Nodejs with MongoDB on server side. In order to make the rendering of polygon layers on the map easy, I would like to create tiles out of GeoJSON object. The output from my Mongo DB for selection of zipcodes looks like this -

[
  {
    properties: { postcode: '47574', locality: 'Goch' },
    geometry: { type: 'Polygon', coordinates: [Array] },
    type: 'Feature',
    _id: 5ed7b221a61a4b2970433932
  },
  {
    properties: { postcode: '52531', locality: 'Übach-Palenberg' },
    geometry: { type: 'Polygon', coordinates: [Array] },
    type: 'Feature',
    _id: 5ed7b221a61a4b2970433933
  },
...
]

I would like to convert this object into GeoJSON from this npm package using const tileIndex = geojsonvt(geoJSON). I am unable to parse this object into valid GeoJSON to create tiles which I can use to render a map on the client-side.


texts: What is the difference between Astar and Dijkstra algorithms in Pgrouting?

texts:
I would like to know what is the difference between A* algorithm and Dijkstra at finding the shortest path at pgrouting.

Which one is faster ,has better performance and returns better results.

Also i would like to know if it is possible for pgrouting to take into account that a road is one way when it computes the shortest path thus not traversing it

Sorry for my english
thanks in advance john


texts: Raster Land Cover Land Use Data added to vector polygon attribute table

texts:
I have a set of 1000 points which I have added a buffer of 100m too. Then I extracted by mask those 1000 points that have the 100m buffer added from a land use land cover (LULC) raster data set. Now I have 1000, 100m circles of raster data. The attribute table for 1000 LULC data still counts the cells of the 1000 different circles as a single item and separates them based on their land class.

I need to have them separated by circle and have the different land classes separated to show how many cells of each land class make up every different circle. I want to be able to see the cells within each circle separately to be able to compare them. I want to add the cell counts of the raster layer to the attribute table of the 100m buffer.


texts: Using Select Layer By Attribute in for loop with ArcPy?

texts:
I'm trying to iterate through a feature class and select each feature one by one by using the OBJECTID field. When I use a single value, the SelectLayerByAttribute works, but how do I compare the OBJEECTID to a variable that increases in count through the loop. The program is then suppose to create a layer for each feature in the feature class. This is what I have so far:

import arcpy
import os
arcpy.env.workspace = "C:\Users\Daimon Nurse\Desktop\DFMPROJECT\DFMPROJECT.gdb"
arcpy.env.overwriteOutput=True
inlayer = "C:\Users\Daimon Nurse\Desktop\DFMPROJECT\DFMPROJECT.gdb\ZonesPolyline"

with arcpy.da.SearchCursor(inlayer, ("OID@", "SHAPE@AREA")) as cursor:
    i = 1
    for row in cursor:
        print("Feature {0} has an area of {1}".format(row[0], row[1]))
        outlayer = "ZonePolygon"
        arcpy.MakeFeatureLayer_management (inlayer, outlayer)
        arcpy.SelectLayerByAttribute_management (outlayer, "NEW_SELECTION", """           "OBJECTID" = 1 """)
        output = r'C:\Users\Daimon Nurse\Desktop\Grounds Project\DFMGROUNDS.gdb'
        outfile = os.path.join (output, i)
        arcpy.CopyFeatures_management(outlayer, outfile)
        print i
        i = i + 1


texts: Mapinfo SQL Query for compare for edited data in tables

texts: I have a table that contain several columns, each has their own data. then, someone have made some change in the data (just a little one data editing). how do i know later what data have been edited/changed using MapInfo SQL Query? thanks.


texts: If statement problem in Field Calculator - Arcgis

texts:
I have a field entitled "Lung_lunc". It contains some "0" values and the rest are values > 0. I want to turn those "0" values, into "0.0001". I have added a new field: "Lung_lunca", and I want to use a IF statement in Field calculator, to do so. I have tried the version in the image, but when I run the script, nothing changes (the values remain the same, as in the previous field).
I have looked over similar questions here, but did not manage to find a solution.

How should I modify the script, in order to work?

P.S. I have made sure that both "Lung_lunc" and "Lung_lunca" are double fields, so they accept decimals.


texts: Moving feature using Python console in QGIS

texts:
I am building a QGIS plugin that can create as much copies of a selected feature as I want and paste it to my required place.I have code to copy and paste feature

select = iface.activeLayer().selectedFeatures()

if self.i==0:
    iface.actionCopyFeatures().trigger()
    self.i=1
else:
    iface.actionPasteFeatures().trigger()

On first click it copies the feature and on second click it paste it and so on. But my problem is I want to paste it on my required place just like I am moving feature from one place to another.

Can you tell me how to move feature and also how to perform this whole task?

I am using this with cursor click so I want to know also how to get feature on my cursor on first click so I can easily paste it on my required place.


texts: Concatenate features in one feature class by the unique ID

texts:
Is there a way to concatenate features in ArcGIS by using ArcPy, or field calculator, to produce the following output?  I am aware of the concatenate tool that Esri created, however it doesn't keep the Features and the Values in the same order.

I need Feature A, to line up with the appropriate Value.  And the tool also condenses multiple values.  If there are two 0's in the output, the tool will get rid of the extra ones.

After talking to Craig.  Here is the python code that I have been testing, this is from the link that he provided in the comments:  I reset the links to my personal folders, however the error that I keep receiving states, "RuntimeError: cannot open" for the link of the new_table.  When I do create a copy of the source table to be used as the new_table, and rename it as table2.dbf, it deletes the information.

import os
import arcpy
arcpy.env.workspace = r"c:\users\Desktop\test"
table = r"c:\users\Desktop\test\GridforPython.dbf" # source_table
new_table = r"table.dbf" # new table name + extension only
tempTable = r"\\in_memory\temp_table"
outPath = r"c:\users\Desktop\test" # path to output location on disk
keyField = "ID_Input"
valField = "Value_Input"

tableDict = {}
arcpy.MakeTableView_management(table, tempTable)
with arcpy.da.SearchCursor(table, [keyField, valField]) as cursor:
    for row in cursor:
        tableDict.setdefault(row[0], []).append(row[1])

insertCursor = arcpy.da.InsertCursor(os.path.join(outPath, new_table), [keyField, valField])
for key, val in tableDict.iteritems():
    insertCursor.insertRow((key, ', '.join(str(item) for item in val)))


texts: Passing data/parameter to Dockable Window in ArcObjects?

texts:
I have a Dockable Window and that dockable window consists of a dataGrid. In order to populate that data grid I want to pass a datatable and some more parameters.

Here is what I did:

private void SetupDockableWindow(DataTable dt)
    {
        if (_dockableWindow == null)
        {
            IDockableWindowManager dockWindowManager = _application as
            IDockableWindowManager;
            if (dockWindowManager != null)
            {
                UID windowID = new UIDClass();
                windowID.Value = _dockableWindowGuid;
                //windowID

                _dockableWindow =
                dockWindowManager.GetDockableWindow(windowID);
                _dockableWindow.Dock(esriDockFlags.esriDockFloat);
            }
        }
    }

How do I pass the required data to it?


texts: Using sieve filter in ERDAS Imagine and/or ArcGIS for Desktop?

texts:
Where can i use a sieve filter for my supervised classification?

Can i do it in Erdas Imagine or in ArcGIS?

I can't find it anywhere.


texts: ST_DWithin use_spheroid argument for data projected in Web Mercator

texts: When calling ST_DWithin in order to fetch data to overlay an area selected on a Web Mercator map, should false be passed to the use_spheroid argument of the call to ST_DWithin, since Web Mercator assumes a perfectly spherical earth? While not geodesically precise, I presume reducing accuracy by using the simplified formulae would actually increase the accuracy of matching the data to the area selected on the map?


texts: Setting source for HTML frame in Print Composer using PyQGIS?

texts:
I have defined an HTML Frame in my composer using:

html_frame = composition.getComposerItemById('html_frame')

How can I set its source using PyQGIS?


texts: What's a good way to implement an 'Attributes to keep' option in a QGIS processing script?

texts:
I'm working on a processing script that selects features by their spatial relationship. The result is a virtual layer:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

# define interface
##Virtual Layers (QGIS 2.x)=group
##Select by spatial relationship=name
##Layer_to_select_from=vector
##Spatial_relationship=selection Intersects;Contains;Within;Touches;Overlaps;Equals
##Intersection_layer=vector
##True_or_false=selection True;False
##Snapping_tolerance=string 1e-5
##Output_layer=string vl_
##Output_unique_identifier_column=string gid
##Output_geometry_column=string geom

from qgis.core import QgsVectorLayer, QgsMapLayerRegistry

# get input parameters from GUI
inlayer_select = processing.getObject(Layer_to_select_from)

spatial_rs = Spatial_relationship
if spatial_rs == 0:
    spatial_rs = 'ST_Intersects'
elif spatial_rs == 1:
    spatial_rs = 'ST_Contains'
elif spatial_rs == 2:
    spatial_rs = 'ST_Within'
elif spatial_rs == 3:
    spatial_rs = 'ST_Touches'
elif spatial_rs == 4:
    spatial_rs = 'ST_Overlaps'
elif spatial_rs == 5:
    spatial_rs = 'ST_Equals'

inlayer_intersection = processing.getObject(Intersection_layer)

t_or_f = True_or_false
snaptolerance = Snapping_tolerance

outlayer = Output_layer
outuic = Output_unique_identifier_column
outgeom = Output_geometry_column

# get geometry type and coordinate reference system from input layer ('Layer to select from')
ingeomtype = inlayer_select.geometryType()
if ingeomtype == 0:
    ingeomtype = 'point'
elif ingeomtype == 1:
    ingeomtype = 'linestring'
elif ingeomtype == 2:
    ingeomtype = 'polygon'

incrs = inlayer_select.crs().authid()
incrs = incrs[5:]

# create query
query_t = """SELECT
    DISTINCT('{0}'.rowid) AS '{4}',
    '{0}'.geometry AS '{5}'/*:{6}:{7}*/
FROM '{0}', '{2}'
WHERE {1}(ST_Snap('{0}'.geometry, '{2}'.geometry, {3}), ST_Snap('{2}'.geometry, '{0}'.geometry, {3}));
""".format( inlayer_select.name(), spatial_rs, inlayer_intersection.name(), snaptolerance, outuic, outgeom, ingeomtype, incrs )

query_f = """SELECT
    '{0}'.rowid AS '{4}',
    '{0}'.geometry AS '{5}'/*:{6}:{7}*/
FROM '{0}' LEFT JOIN '{2}'
ON {1}(ST_Snap('{0}'.geometry, '{2}'.geometry, {3}), ST_Snap('{2}'.geometry, '{0}'.geometry, {3}))
WHERE '{2}'.rowid IS NULL;
""".format( inlayer_select.name(), spatial_rs, inlayer_intersection.name(), snaptolerance, outuic, outgeom, ingeomtype, incrs )

if t_or_f == 0:
    query = query_t
elif t_or_f == 1:
    query = query_f

# run query and add layer to map
vlayer = QgsVectorLayer( "?query={}".format(query), outlayer, "virtual" )
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

Here is a short screencast if someone wants to see the script in action: http://hub.qgis.org/attachments/10767/select_by_spatial_relationship.mp4

At the moment the GUI misses an 'Attributes to keep' option. As far as I know there is no parameter type to list the attributes of an input layer as checkboxes (see http://docs.qgis.org/2.14/en/docs/user_manual/processing/scripts.html and http://hub.qgis.org/issues/16348). Would an input field (type string), that is designed to accept a comma separted list, do the job? Can anyone help me to implement an 'Attributes to keep' option in a smooth way?


texts: Use PHP to convert Texas State Plane to decimal

texts:
I'm trying to find a PHP function to convert a state plane coordinate to decimal.

For example,

NAD_1983_stateplane_texas_south_fips_4205_feet (EPSG:102741)

to

WGS84 (EPSG:4326)... or NAD83


texts:
I'm trying to get all the nodes of one street in a district in Berlin.

I've written the following expression, but it returns all Holteistraße in Germany. Anyone know how I would write a correct syntax?

[out:json];

(
area["Berlin, Friedrichshain"];
way[highway][name="Holteistraße"];
);

out body;
>;
out skel qt;


texts: I am starting to develop a plugin for Qgis using python, which i am very new at, and would like to view some results of variables by displaying a message in QGIS. How can I do this?


texts:
I'm trying to write a script to copy multiple feature classes from two different databases into a single database, but I want to be able to specity what features I actually want. I know I need to use da.Walk to get through the first files but am stumped on how I can actually do this. Any thoughts?

here is what i have so far.

`import arcpy, os, string

workspace= ""

for dirpath, dirnames, filenames in arcpy.da.Walk(workspace, datatype="FeatureClass"):
    for filename in filenames:`


texts: Assessing impact of spatial scale on derivative calculations and associated kernel size?

texts:
Full question: What is the most appropriate way to assess the impact of spatial scale on derivative calculations associated with typical 3x3 moving window kernel operations?

This relates to a couple of existing questions:

Finding local maxima with variable window search using ArcGIS Desktop? and
Using variable window ("kernel") size for raster processing in QGIS?

but the problem is not the same.

The issue presented here relates to how to best account for spatial variability whilst treating element positions within a kernel consistently whilst also minimising potential uncertainty inherent of some possible processing steps.

Background

For many terrain analysis algorithms, a moving 3x3 kernel is often implemented of the form:

[i-1, j-1] [i,j-1] [i+1,j-1]
[i-1,j]     [i,j]  [i+1,j]
[i-1,j+1] [i,,j+1] [i+1,,j+1]

This works for many implementations such as slope calculations where each pixel is equally weighted. The jump I'm not able to make is what happens where we are interested in implementing a function over a larger area e.g. to account for spatial dependence of the result?

As far as I can see, there are two ways of addressing such a problem:

Increase the kernel pixel dimensions whilst leaving the underlying raster at its original cell resolution
Keeping the kernel pixel dimensions the same but resampling the raster so that the kernel is by default covering a larger area

Option 1 - increase kernel dimension

Where all values within the kernel are operated on equally e.g. a mean is taken, then the kernel can simply be expanded so 3x3 >> 5x5 >> 7x7 etc. However where specific locations have specific importance e.g. when calculation Horn's (1981) finite difference estimator (east-west gradient):

(Horn's (1981) finite difference estimator taken from Burrough and McDonnell, (1998))

bigger kernels can introduce hurdles.... For example, it would be easy to over-weight values in the corners of the kernel as opposed to the horizontal/perpendicular locations. Thus, how you specifically weight pixels is something to consider.

Option 2 - resample raster using fixed kernel pixel dimension

Resampling the raster will force the window to cover a larger area (a pixel is still a pixel after all) however by resampling our underlying raster we:

introduce uncertainty associated with the chosen resampling technique
Lose some of the higher frequency components (perhaps desirable....)
Remove the ability to cleanly assign any derived values from the moving window to a raster of the same dimensions of the original

So....

Methods do exist considering multi-scale changes e.g. Topographic Position Index (see poster here) but this doesn't address my conundrum.

I've played with spatial scale in the past with various image processing applications but these have always taken a mean value from a kernel of varying size - they have not been concerned with the specific locations of pixels within the kernel.

I am certainly not the first to think about this problem - does anyone have any references they can point to that deal with the issues outlined above?


texts: How to obtain a WKT from a DotSpatial IFeature geometry?

texts:
I am using DotSpatial 1.9 and I have a feature (sourced from a shapefile).

According to the docs, I should be able to call: ToString() on IFeature.

    //
    // Summary:
    //     Returns the Well-known Text representation of this Geometry. For a definition
    //     of the Well-known Text format, see the OpenGIS Simple Features Specification.
    //
    // Returns:
    //     The Well-known Text representation of this Geometry.

But what I getting is the name of the type:

feature.ToString()
"DotSpatial.Data.Feature"

What am I missing?


texts: Problems with PostGIS St_closestpoint function

texts:
I'm trying to find the closest point in a line from an insertion point closest to it. To do it, in the middle of the function I use
select id from rede_viaria.t_rviaria_via vias
where
  st_intersects(
    st_closestpoint(
      vias.the_geom,
      (st_transform(
        ST_GeomFromEWKT(
          'SRID=3857;POINT(-922235.81565421 5053077.824914)'
        ), 3763)
      )
    ),
    vias.the_geom
  );

The problem is that I'm receiving 1200 rows instead 1 based on the intersection. Could anyone please help? What I'm doing wrong?