texts: automatically georeference QgsGeometry.fromPoint features on creation

texts:
Is there a way to automatically georeference fromPoint features upon creation?

I use this layer URI:

      self.pinLayer =  QgsVectorLayer(
          "Point?crs=epsg:3857&field=id:integer&field=description:string(120)&field=x:double&field=y:double&field=pluscode:string(20)&index=yes",
          "Pins",
          "memory")

I then set the pins using the QgsGeometry.fromPoint method:

        feature = QgsFeature()
        feature.setGeometry(QgsGeometry.fromPoint(point))

Which works fine as you can see:

But the points are not georeferenced.
Is there a way to add / do that right during creation?

I've found this old entry but its for raster layers and it also requires some user interaction (setting ground control points)
I want to do it programmatically


texts: tileLayers grouped in layerGroup not displaying

texts:
I am trying to combine 3 tileLayers (Esri map services for imagery, boundaries, roads) into one layerGroup so that it can be chosen via the layers control.

This jsfiddle shows the three tileLayers initially loaded individually ok.  But the layerGroup doesn't work (the map is gray) when 'Grouped' is chosen via the L.control.layers.


texts: What is Oracle/ArcSDE doing to my simple WKT Conversions?


texts:
I am working with rule based labels on a line layer and need to move the labels along different axes along the line. I want to move the label 2 along the line a few millimeters.
For the label 2, the placement is a curve, and the position is above the line. The labels 1 and 2 will repeat after 15 millimeters, and I need the label 2 to be offset from the label 1 along the line by a few millimeters. There is a distance parameter, but it is responsible for the offset not along the line, but across.

Сurrent setting for the first label
Current setting for the second label


texts: Add multiple geotiff files in a geoserver coverage store

texts:
I'm using the Geoserver REST interface to create new coverage stores and add geotiff files to it.

I have a zip file with multiple tiff files in it, one for each band, and would like to have them all stored under the same coverage store. The main goal is then to be able to combine bands using the coverage view.

My problem: the coverage store seems to only accept one tiff file at a time.

How can I add multiple tif files to a coverage store?


texts: How to merge elevation from a raster image onto a shape road network with QGIS?

texts:
I have a shape file with Stockholms road network and a raster image with elevation data from the same area.

Is there any way to merge the two somehow, either by adding the slope direction as positive or negative numeric values through an extra column in the attribute data or in any other magic way?

The purpose would be to calculated shortest path with regard to the elevation, to be able to avoid unnecessary uphill climbs.

I work with QGIS and pgRouting so an answer compatible with that would preferable, but other solutions are welcome as well!


texts: I have groups of vector points (GPS points) in a number of layers in QGIS.  I would like to connect the points in each layer into polygons in each individual layer.  I have look hard, but can not find away within QGIS to do this ?


texts:
I have a ArcObjects 10.2 vb.net Custom Add-in that is a dockable window.  The dockable window displayes a windows form that has listboxes, textboxes, panels, tab panels, etc. that the user interacts with.   I'm trying to get it so that when the window is not docked it autosizes like a windows form.  I want certain controls on the form to grow/shrink when the dockable window is resized.  I'm struggling to get this to work and wonder if it is even possible.  If I anchor the controls on all sides they overlap each other in the dockable window. I've tried multiple combinations of anchors, autosizing modes on both the dockable window and the form.

I have looked at resizing the window here:
http://help.arcgis.com/en/sdk/10.0/arcobjects_net/componenthelp/index.html#/Resize_Dockable_Window_Snippet/004900000077000000/

Moving the window here:
http://help.arcgis.com/en/sdk/10.0/arcobjects_net/componenthelp/index.html#/Report_Dockable_Window_Position_Snippet/0049000000qq000000/

And creating a fixed size here:
https://gis.stackexchange.com/questions/143623/creation-of-fixed-size-dockable-window

I haven't found anything that shows how this can be accomplished.  Is there away to get the dockable window to act like a windows form with the same autosizing properties and still dock to ArcMap?


texts: Open data cube GeoTIFF indexing from MinIO/S3

texts:
I have large amounts of existing GeoTIFF data and I would need to find example how to index those to open data cube. Currently those reside in MinIO (S3 compliant store). In my understanding I would need to convert these to COG (that I can probably manage) and then I would need to generate STAC file to ingest those. Any pointer to tutorial or examples?

Thanks for the tips. However I think s3-to-dc won't work with Minio as there is no way to pass endpoint.


texts: Is there any link form where I can download Sentinel 2 Level 2A products for Asian Countries or it has to be done through Sentinel Toolbox only?


texts:
I made a shapefile from a DEM raster of Ethiopia, and I want to find the total area of the shapefile. I have changed the CRS of all files to Africa_Albers_Equal_Area_Conic (ESPG: 102022), and the units are meters. I have tried using the field calculator, but the answers are very small. I have also used the identify features tool, and the answers are the same as the field calculator's.

I also used "count" in the field calculator to get the number of pixels, but these sometimes seem very small too. The pixel size is 0.000277778,-0.000277778, but I don't know what the units are, so I am not sure how to use that to find the area.

I tried using the measure area tool where you draw a polygon, but this gave a very small answer too.

I read there might be a way to do it with GRASS, but I have a problem with my GRASS plugin which I don't understand.


texts:
Is there a good tool / technique to display ArcSDE-metadata in a webbrowser?

We want to provide user without GIS-knowledge with a user friendly tool to explore our SDE metadata catalog. We are using the inspire-profile for metadata input. They just need to explore the properties of the metadata, they don't need to see the data itself.

I am currently experimenting with geoportal to allow user without arcgis explore our metadata catalog. But i feel like it is to complex for non-gis users, with all the categories etc.

I have also looked at geonetwork, but there seems to be no way to automate the update from database to  server and i also had issue's with the import of metadata xml from exported form arccatalog.


texts: Convert GeoJSON to GML using an old GeoTools library

texts:
I'm trying to convert a GeoJSON file to GML 3.1.1 using an old org.geotools lib. geotools lib are 2.7.0 version. Sadly, upgrading is not an option. Here is the relevant part of my code:

FeatureCollection geoJsonDocument = featureJSON.readFeatureCollection(myGeoJsonInputStream));
...
Configuration config = new org.geotools.gml3.GMLConfiguration();
Encoder encoder = new Encoder(config);
encoder.getNamespaces().declarePrefix("myorgns", "http://myorg.org");
encoder.setEncoding(StandardCharsets.US_ASCII);
encoder.encode(geoJsonDocument, new QName("http://www.opengis.net/gml", "FeatureCollection"), out);

If I print the content of out variable I find that GeoJSON elements in GML all have "null" prefix. Here is a snippet of the content:

<gml:FeatureCollection xmlns:sch="http://www.ascc.net/xml/schematron" xmlns:gml="http://www.opengis.net/gml" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:myorgns="http://myorg.org">
<gml:featureMembers>
    <null:feature gml:id="feature-0">
        <gml:boundedBy>
            <gml:Envelope>
                <gml:lowerCorner>3.4686311 51.1304626</gml:lowerCorner>
                <gml:upperCorner>3.4686311 51.1304626</gml:upperCorner>
            </gml:Envelope>
        </gml:boundedBy>
        <null:pk_track>729</null:pk_track>
        <null:time_ISO8601>2017-01-10T05:13:06+01:00</null:time_ISO8601>
        <null:time_epoch>1484021586000</null:time_epoch>
...

When I try do encode the geoJsonDocument I got an exception because encoded element prefixes are "null":

Exception in thread "Thread-30" java.lang.RuntimeException: org.xml.sax.SAXParseException; lineNumber: 2; columnNumber: 238; The prefix "null" for element "null:feature" is not bound.
at org.geotools.xml.StreamingParser$1.run(StreamingParser.java:248)
at java.lang.Thread.run(Thread.java:748)
Caused by: org.xml.sax.SAXParseException; lineNumber: 2; columnNumber: 238; The prefix "null" for element "null:feature" is not bound.
at org.apache.xerces.util.ErrorHandlerWrapper.createSAXParseException(Unknown Source)
at org.apache.xerces.util.ErrorHandlerWrapper.fatalError(Unknown Source)
at org.apache.xerces.impl.XMLErrorReporter.reportError(Unknown Source)
at org.apache.xerces.impl.XMLErrorReporter.reportError(Unknown Source)
at org.apache.xerces.impl.XMLErrorReporter.reportError(Unknown Source)
at org.apache.xerces.impl.XMLNSDocumentScannerImpl.scanStartElement(Unknown Source)
at org.apache.xerces.impl.XMLDocumentFragmentScannerImpl$FragmentContentDispatcher.dispatch(Unknown Source)
at org.apache.xerces.impl.XMLDocumentFragmentScannerImpl.scanDocument(Unknown Source)
at org.apache.xerces.parsers.XML11Configuration.parse(Unknown Source)
at org.apache.xerces.parsers.XML11Configuration.parse(Unknown Source)
at org.apache.xerces.parsers.XMLParser.parse(Unknown Source)
at org.apache.xerces.parsers.AbstractSAXParser.parse(Unknown Source)
at org.apache.xerces.jaxp.SAXParserImpl$JAXPSAXParser.parse(Unknown Source)
at org.apache.xerces.jaxp.SAXParserImpl.parse(Unknown Source)
at javax.xml.parsers.SAXParser.parse(SAXParser.java:195)
at org.geotools.xml.StreamingParser$1.run(StreamingParser.java:244)
... 1 more

How can I solve this?


texts: Modify expression so it can act in a group of points with a value in common

texts:
Can somebody help me to modify this expression so it can act in a group of points with a value in common? As it is right now it numbers the points extracted from a line based on the northernmost point clockwise, but it takes all the points from the layer. I wished it did something as in the image bellow.
if (
    vertex_index <
    aggregate(
        @layer,
        'max',
        vertex_index,
        $y = maximum( $y)
    ),
    vertex_index + ( aggregate (
        @layer,
        'count',
        vertex_index
    ) - aggregate (
            @layer,
            'max',
            vertex_index,
            $y = maximum( $y)
    )),
    vertex_index - aggregate (
        @layer,
        'max',
        vertex_index,
        $y = maximum( $y)
    )
)


texts: Manage files in geodatabase for ArcPy loop?

texts:
I have a piece of code consisting of hundreds of lines that contain multiple processes (snapping, feature to polygon, erase, etc.) and I am going to make a for loop that will take care of many files that can be found in a geodatabase.

How to keep track of all the files and their copies in the cases where you need to maintain copies in operations like Snap_edit where it modifies the original and use only the files needed for the processes that will follow?

It looks chaotic and would like to hear ideas about how to manage the files (original files, copies, outputs of mid-processes etc.) in the geodatabase so they can be used for the processes in a for loop.


texts: I'm looking for a brief history of a 256x256 web map tile, in particular who used it first (came up with the idea) and who popularized it.


texts:
I have a Windows Vista (32 bit) computer.  I tried installing and reinstalling (by uninstalling the software first) QGIS 1.8 a few times then went to version 1.7.4 and still encountering the same problem.

I don't know what the problem is since I am downloading the Windows version of the software.  Could anyone help me with this issue?


texts: Does anybody know of a software that can project VIIRS EDR products. I used ENVI 5.3 to project vegetation index products available on NOAA CLASS but ENVI project the images in a weird way. The two sides of the swath get projected fine but the middle of swath is missing in the projected image. Please let me know if you know any more efficient way to project these data.


texts:
I currently use Quantum GIS 1.7.4 at work and I have trouble to open Google layers extension. In fact, when I start QGIS 1.7.4 the following message appears and I cannot use Google layers extension. Would really appreciate if anyone can help me sort this problem out.

Traceback (most recent call last):

  File "C:/PROGRA~2/Quantum GIS Wroclaw/apps/qgis/./python\qgis\utils.py", line 147, in startPlugin
    plugins[packageName].initGui()
  File "C:/Users/JICA/.qgis/python/plugins\GoogleLayers\GoogleLayers.py", line 33, in initGui
    self.sessionid = gmaps.get_sessionid()
  File "C:/Users/JICA/.qgis/python/plugins\GoogleLayers\gmaps.py", line 21, in get_sessionid
    sessionid = t.strftime("%d%m%Y_%H%M%S%s")
ValueError: Invalid format string

Version de Python :
    2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]

Version de QGIS :
    1.7.4-Wroclaw Wroclaw, 411aff6

Chemin vers Python :

['C:\\Users\\JICA\\.qgis\\python\\plugins\\GeoCoding', 'C:\\Users\\JICA\\.qgis\\python\\plugins\\elevation', 'C:/PROGRA~2/Quantum GIS Wroclaw/apps/qgis/./python', 'C:/Users/JICA/.qgis/python', 'C:/Users/JICA/.qgis/python/plugins', 'C:/PROGRA~2/Quantum GIS Wroclaw/apps/qgis/./python/plugins', 'C:\\PROGRA~2\\Quantum GIS Wroclaw\\bin\\python27.zip', 'C:\\PROGRA~2\\Quantum GIS Wroclaw\\apps\\Python27\\DLLs', 'C:\\PROGRA~2\\Quantum GIS Wroclaw\\apps\\Python27\\lib', 'C:\\PROGRA~2\\Quantum GIS Wroclaw\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~2\\Quantum GIS Wroclaw\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~2\\Quantum GIS Wroclaw\\apps\\qgis\\bin', 'C:\\PROGRA~2\\Quantum GIS Wroclaw\\apps\\Python27', 'C:\\PROGRA~2\\Quantum GIS Wroclaw\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~2\\Quantum GIS Wroclaw\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~2\\Quantum GIS Wroclaw\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~2\\Quantum GIS Wroclaw\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~2\\Quantum GIS Wroclaw\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~2\\Quantum GIS Wroclaw\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\Users\\JICA\\.qgis\\python\\plugins\\fTools\\tools']

Thanks for your help.


texts: QGIS processing algortihm crashes in python console but works manualy

texts:
I am loading two layers from PostGIS database and I am trying to run 'qgis:intersection' algorithm with them. Execution of the algorithm through the python console initially works but crash after 2-3%. Loading the layers through the console and then manually starting the algorithm works all the way though. I'm guessing that there is further processing of the layer that needs to be before putting them through the intersection algorithm. Here is an example of my script.

import processing

sql1 = "select system, gid, surface, st_collectionextract(st_makevalid(geom),3)::geometry(MultiPolygon, 3047) as geom from dobrich_ias_final"

uri1 = QgsDataSourceURI()
uri1.setConnection("****", "5435", "****", "****", "****")
uri1.setDataSource("" ,'(%s)' % sql1 , "geom","","gid")

ias = QgsVectorLayer(uri1.uri(),"ias","postgres",True)

if not ias.isValid():
    print "Layer %s did not load!" % ias.name()

sql2 = "select sub_id, gid, pop,tot_area, area_1, area_2, area_3, st_collectionextract(st_makevalid(geom),3)::geometry(MultiPolygon, 3047) as geom from dobrich_dobrich_s_subcatchments_result_m2_30"

uri2 = QgsDataSourceURI()
uri2.setConnection("****", "5435", "****", "*****", "*****")
uri2.setDataSource("", '(%s)' % sql2 , "geom","", "gid")

subs = QgsVectorLayer(uri2.uri(),"subs","postgres",True)

if not subs.isValid():
    print "Layer %s did not load!" % subs.name()

QgsMapLayerRegistry.instance().addMapLayers([ias,subs],False)

i_res=processing.runalg('qgis:intersection',ias,subs,False,None)


texts:
I'm downloading more than 700 tiled SHP from this site:
"http://www502.regione.toscana.it/geoscopio/cartoteca.html#"
It is possible to download automatically all tiles from QGIS or ARCGIS?

The site provide also a qgis an interactive qgis project, like this:

This is wms link:
http://www502.regione.toscana.it/wmsraster/com.rt.wms.RTmap/wms?map=wmscartoteca&map_resolution=91&


texts: QGIS Deployment in Citrix XenApp - How to Remove Pyramids from the image "properties"

texts:
I'm working on deploying QGIS 2.4 into a XenApp 6.5 desktop environment.

We've tested QGIS to see if it will run in our environment on a test server and have found that building Pyramids with aerial photos for a map causes the user's session to grab as much memory as it wants. is there either a way to limit the memory a session can grab, or disable the building Pyramid feature?

We load balancing our XenApp environment, we have 2 different memory configurations running in this farm, we have 48GB and 32GB servers. If we ended up with 2 QGIS users on 1 32GB server and 1 of them ran the Build Pyramid command, it could potentially cause the server to fall over.

Apologies if this is the wrong place to post.  Please point me in the right direction if that's the case.


texts: How to read slope from HYDRO1k database for Asia

texts:
The Hydro1k database for Asia is provided by USGS. It is available here as GTOPO30 HYDRO 1K.

I tried to read the slope values for the Kumamoto area in Japan but the results all came out as 55537(same as for the sea). The slope is shown different when I open it in QGIS.

Here is the prj file:

PROJCS["US National Atlas Equal Area",GEOGCS["Unspecified datum based upon the Clarke 1866 Authalic Sphere",DATUM["D_Sphere_Clarke_1866_Authalic",SPHEROID["Clarke_1866_Authalic_Sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_origin",45],PARAMETER["central_meridian",100],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]

Here is my code:

import osgeo
from osgeo import gdal,osr
import ogr
geo = gdal.Open('/home/Downloads/gt30h1kas/as_slope.bil') #slope file

# Read the projection from prj file
ds = gdal.Open('/home/Downloads/gt30h1kas/as_dem.prj')
prj=ds.GetProjection()
geo_ref=osr.SpatialReference(wkt=prj)

# Reference system for input coordinates
point_ref=ogr.osr.SpatialReference()
point_ref.ImportFromEPSG(4326)

ctran2=ogr.osr.CoordinateTransformation(point_ref, geo_ref)

arr = geo.ReadAsArray()

geotransform = geo.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]

# Kumamoto (32.8031° N, 130.7079° E)
[lon,lat,z]=ctran2.TransformPoint(130.7079, 32.8031)
print(lon,lat)
offx = int( (lon-originX)/pixelWidth )
offy = int( (lat-originY)/pixelHeight )
print(offx, offy, arr[offx][offy])

Here is the output:

2810038.798891878 -837540.6774932265
7165 7281 55537

Please help me find my mistake.


texts: How to access the polygon vertices in ArcObjects

texts:
I am able to retrieve a list of polygons from the opened map layer in ArcGIS using the following C# code:

        IMxDocument pMxDoc;
        pMxDoc = (IMxDocument)ArcMap.Application.Document;
        IMap pMap;
        pMap = pMxDoc.FocusMap;

        IFeatureLayer2 pFLayer;
        pFLayer = (IFeatureLayer2)pMap.get_Layer(0);

        // Retreve the feature class.
        IFeatureClass featureClass = pFLayer.FeatureClass;

        // Retrieve the number of polygons in the map.
        IQueryFilter queryFilter = new QueryFilterClass();
        int numberOfPolygons = featureClass.FeatureCount(queryFilter);

        for (int i = 0; i < numberOfPolygons; i++)
        {
            IFeature feature = featureClass.GetFeature(i);

            // Get polygon.
            IGeometry geometry = feature.Shape;
            IPolygon polygon = (IPolygon)geometry;
        }

How to access the vertices (coordinates) of the rings of those polygons?


texts: Determining homogenity of flatroofs based on digital surface model / slope layer using QGIS

texts:
For the city of Berlin I have a digital elevation model as a 1x1 raster dataset from what I counted the slopes and I have shapefile data containing the outlines of buildings. My goal is to find flatroofs with preferably little fractionation.
Therefore I eliminated all pixel values with a slope higher than 3 degrees. An example of the result is seen in the picture. Now I would like to assess the "quality" of the roofs, so where the most area of flat pixels is "connected" like it is the case with the roof of the right building in the example vs. the more destructed roof on the left. As a result I would like to classify the roofs in little/medium/much destruction.
How do I approach this using QGIS?


texts: Using Orfeo Toolbox object based classification in QGIS tutorial/workflow?

texts:
I have a set of aerial photos of vegetation which needs to be classified. Instead of using manual interpretation, I want to try using Object-Based classification based on the Orfeo Toolbox in QGIS Processing.

However, I cannot find any tutorial that explains how to do it in QGIS. In the Orfeo guide, the steps are

Image segmentation (the whole or only parts of it);
Image to LabelObjectMap (a kind of std::map) transformation;
Eventual relabeling;
Attribute computation for the regions using the image before segmentation:
Object filtering
LabelObjectMap to image transformation.

But I cannot find any specific step-by-step, which geo algorithm to use.
I tried to follow LSMC (Large-Scale Mean-shift Classification) workflow, and also run Segmentation (Meanshift). The LSMC result was quite good, but I have no idea what to do next to label the object (or make a training object) and classify it.

I tried the TrainImageClassifier(SVM), TrainORGLayerClassifier, ORGLayerClassifier, but there were no outputs came out. I tried to find some tutorials, but could not find any. I am totally confused about the workflow.


texts:
I am trying to migrate a plugin to QGIS3 from QGIS2 and I am having issues with an API.

The original python has QgsMessageBar.MessageLevel but https://qgis.org/api/api_break.html says it has been moved to Qgis::MessageLevel. I tried replacing it in the code and got an invalid syntax and I couldn't find a MessageLevel example in py3.

What does it mean when something has been moved to Qgis::...? How do I find the proper syntax and use?


texts: CartoDB restricts uploading only up to 250 columns

texts: I had developed an application using cartodb.js which transforms a csv file uploaded in cartodb database to geometries in map. The application works fine when there are 230 columns in csv file and the upload is not success if i tried to import 309 columns into cartodb. Is there any other way to increase the column count or any alternate solutions?


texts: Arcgis 3D map development when display PictureMarkerSymbol in sceneview

texts:
When I am displaying PictureMarkerSymbol it gives below exception in console.

> DOMException: Failed to execute 'texImage2D' on
> 'WebGLRenderingContext': The cross-origin image at "image url here"
> may not be loaded.(…) "Error: Failed to execute 'texImage2D' on
> 'WebGLRenderingContext': The cross-origin image at
> https://webapps-cdn.esri.com/Apps/MegaMenu/img/logo.jpg may not be
> loaded.
>     at Error (native)
>     at Object.makePotTexture (http://js.arcgis.com/4.0/esri/views/SceneView.js:244:181)
>     at Object.texImage2D (http://js.arcgis.com/4.0/esri/views/SceneView.js:241:437)
>     at b.loadGl (http://js.arcgis.com/4.0/esri/views/SceneView.js:238:366)
>     at e.aquire (http://js.arcgis.com/4.0/esri/views/SceneView.js:269:453)
>     at f.aquireIfNotUndefined (http://js.arcgis.com/4.0/esri/views/SceneView.js:84:404)
>     at Object.f.singleTextureGLMaterialConstructor (http://js.arcgis.com/4.0/esri/views/SceneView.js:86:1)
>     at new x (http://js.arcgis.com/4.0/esri/views/SceneView.js:667:213)
>     at h.aquireExt (http://js.arcgis.com/4.0/esri/views/SceneView.js:418:406)
>     at h.aquire (http://js.arcgis.com/4.0/esri/views/SceneView.js:417:465)
>     ----------------------------------------
>     rejected at c (https://js.arcgis.com/4.0/init.js:154:313)
>     at e (https://js.arcgis.com/4.0/init.js:154:65)
>     at f (https://js.arcgis.com/4.0/init.js:153:349)
>     at b.resolve (https://js.arcgis.com/4.0/init.js:155:406)
>     at _loadedHandler (http://js.arcgis.com/4.0/esri/views/SceneView.js:1337:242)
>     at http://js.arcgis.com/4.0/esri/views/SceneView.js:793:83
>     at q.done (http://js.arcgis.com/4.0/esri/views/SceneView.js:793:298)
>     at _doneLoadingCallback (http://js.arcgis.com/4.0/esri/views/SceneView.js:789:529)
>     at z (http://js.arcgis.com/4.0/esri/views/SceneView.js:798:307)
>     at HTMLImageElement.d.onload (http://js.arcgis.com/4.0/esri/views/SceneView.js:790:332)
>     ---------------------------------------- Error
>     at then.c.then (https://js.arcgis.com/4.0/init.js:156:177)
>     at https://js.arcgis.com/4.0/init.js:177:128
>     at _prepareImageResources (http://js.arcgis.com/4.0/esri/views/SceneView.js:1122:284)
>     at _prepareResources (http://js.arcgis.com/4.0/esri/views/SceneView.js:1119:301)
>     at constructor (http://js.arcgis.com/4.0/esri/views/SceneView.js:1133:34)
>     at new <anonymous> (https://js.arcgis.com/4.0/init.js:45:148)
>     at Object.make (http://js.arcgis.com/4.0/esri/views/SceneView.js:1116:186)
>     at constructor (http://js.arcgis.com/4.0/esri/views/SceneView.js:1113:287)
>     at new <anonymous> (https://js.arcgis.com/4.0/init.js:45:148)
>     at e.getOrCreateGraphics3DSymbol (http://js.arcgis.com/4.0/esri/views/SceneView.js:1069:390)"


texts: Finding the outer ends of lines in a line feature class

texts: I have a polyline feature class and need to turn the ends of the lines into points. However, the lines are made up of segments that are not quite contiguous, with slight gaps between many of the segments. In effect, the lines are dashed, with varying lengths of the gaps and segments. Is there any way, in ArcGIS, to find only the points that are at the ends of the dashed lines, not at the starts and ends of the "dashes"?


texts: Building Raster Attribute Table (RAT) from GeoPandas dataframe with rasterio or GDAL

texts:
I am trying to figure out how to generate a RAT for a 1-band raster with rasterio or GDAL given a GeoPandas dataframe.

My GeoPandas dataframe has 4 columns (X, Y, PERC_VALUE, LANDUSE_LABEL) which I would need to use to build my raster as so:

X, Y defines the coordinates of the pixel centroid
PERC_VALUE is the value which has to be stored in the pixels
LANDUSE_LABEL has information associated to the pixels which I need to store in the RAT.

The Pandas dataframe that I am using to build the GeoPandas dataframe looks like this:

landuse value   Y   X
FOREST  0.0005  1386050     4559950
CITY    0.0001  1386250     4559650
FOREST  0.0009  1386250     4559750
FOREST  0.0025  1386250     4559850
CITY    0.0059  1386250     4561950
...

I am able to create a raster to store the coverage percentage values at each corresponding X,Y pixel (code is shown at the bottom of the question). What I don't know is how can I associate each pixel to the corresponding LANDUSE_LABEL information in a RAT?

The idea is that when I open the final result with (e.g.) QGIS (if it can read RATs), I would be able to symbolize the raster per landuse label, meanwhile being able to access the percentage values of the pixels.

Is this possible with rasterio or GDAL and if so how? Do I have to use a specific raster extension for this aim (I am currently using GeoTIFF)?

Code I am using to write the PERC_VALUE with rasterio as GeoTIFF:

# convert Pandas df to GeoPandas DF (Points)
gdf = GeoDataFrame(
    dataframe.drop(['X', 'Y'], axis=1),
    crs={'init': 'epsg:3035'},
    geometry=[Point(xy) for xy in zip(dataframe.X.astype(int), dataframe.Y.astype(int))])

# write the raster
with rasterio.open(r'output.tiff',
               'w+',
               driver= 'GTiff',
               height= arr.shape[0],
               width= arr.shape[1],
               count= 1,
               dtype= numpy.dtype('float32'),
               crs= SR,
               transform= transform) as out:

    out_arr = out.read(1)

    # this is where we create a generator of geom, value pairs to use in rasterizing
    shapes = ((geom,value) for geom, value in zip(gdf['geometry'], gdf['value']))
    burned = features.rasterize(shapes=shapes, fill=0, out=out_arr, transform=transform)
    out.write(burned, 1)

    out.close()


texts:
I'm downloaded OpenLayers v6.3.1 libraries and successfully displayed a simple map on my local server, and nom I'm trying to load a GeoJSON file as a new vector layer into this map.

But for some reason it's not loading...

Here is what it shows in Chrome :

And here is the code.

<!doctype html>

    <link rel="stylesheet" href="css_cyclo_temp.css">
    <link rel="stylesheet" href="libraries/ol.css">
    <title>TestMap</title>

I'm quite on beginner on this OpenLayer thing... Thanks a lot ^^'


texts:
I am tying to use pgrouting in order to calculate the distance that can be passed within a certain time period. Someone has provided me an updated OSM dataset with speeds. From there I am trying to use pgrouting in order to calculate from a certain point on that net to every line segment the distance.

Pgrouting is documented, however there are not that many examples out there. While I think it is an absolute amazing tool. Hence my question here.

------- UPDATE

select pgr_createTopology('temp.bikenet_tilburg_4326', 0.000001, 'geom', 'gid');

= ok

SELECT pgr_analyzegraph('temp.bikenet_tilburg_4326', 0.000001, 'geom', 'gid');

= ok

SELECT * FROM pgr_drivingDistance('SELECT gid as id, source, target, time as cost FROM schema.table', 1::bigint, 12::bigint, false)

The drivingdistance only shows how many line segments i had to pass, which is the costs. But I want it to count the time each segment takes. Think I can built that on top of it but there might be a better option inside.


texts: Creating field with specified slope projecting from 3D polygon?

texts:
I am trying to figure out how to create a elevation raster from a known location (3D polygon) that projects at a specified slope (1%) for a specified distance(1-2 km).

Does anyone have any idea how to do this?


texts: OpenLayers Google v3 and other layers drift apart

texts:
I'm stuck. My code works with OSM, but I need it to work with Google as well.  The problem is simply that Google v3 and my other layers are not playing well together at all.  They don't zoom together. My problem is specifically with OL 2.13 and 2.12, and it didn't start happening till about a month ago.

In the case of Google satellite, the layer doesn't zoom at all unless I pan the map.  In the case of the others (hybrid, terrain, roads), they zoom, but the google layer and the other layers on the map zoom up to half a second apart from each other, and the layer "drifts" between zoom levels somewhat.  It shows up markedly off by a few hundred feet and gets worse as I zoom in or out from the initial zoom level. (It appears to be correct at the initial zoom level).

The code that reproduces this is simple and derived directly from the WMS and Google v3 examples, but they appear not to work together well. Vector layers have the same behavior:

$(function() {
    var in_options = {
        internalProjection: new OpenLayers.Projection('EPSG:900913'),
        externalProjection: new OpenLayers.Projection("EPSG:4326")
    };

    var wkt = new OpenLayers.Format.WKT(in_options);
    var feature = wkt.read("POLYGON ((-78.9950752916525687 35.5198238058974169, -78.9950752916525687 36.0750238562958074, -78.2565296570171398 36.0750238562958074, -78.2565296570171398 35.5198238058974169, -78.9950752916525687 35.5198238058974169))");
    var bounds = feature.geometry.getBounds();

    var wms = new OpenLayers.Layer.WMS('Zoning', "/ga_resources/wms/", {

       styles: "styles/new-style",
       layers: "catalog/states/nc/counties/wake/zoning",
       transparent: true
    });

    var map = new OpenLayers.Map({ div:'map', projection: new
        OpenLayers.Projection("EPSG:900913")
    });

    var layers = [
      new OpenLayers.Layer.Google("Google Physical", {type: google.maps.MapTypeId.TERRAIN}),
      new OpenLayers.Layer.Google("Google Streets",{numZoomLevels: 20}),
      new OpenLayers.Layer.Google("Google Hybrid",{type: google.maps.MapTypeId.HYBRID, numZoomLevels: 20}),
      new OpenLayers.Layer.Google("Google Satellite",{type: google.maps.MapTypeId.SATELLITE, numZoomLevels: 22}),
      wms
   ];

   map.addLayers(layers);
   map.addControl(new OpenLayers.Control.LayerSwitcher());
   map.zoomToExtent(bounds);
});


texts: Pass information in external command line software from a python plugin file

texts:
it would be quite difficult to explain, I'll do my best..

I'm writing a plugin for QGIS (with the plugin builder, so the folder structure is straightforward). This plugin calls other software (that of course have to be installed on the computer) that do some mathematical models.

Now, the problem is that I'm able to call these software (all these software work from the command line) from the python plugin but at this point I have to pass to the opened command line, the commands required by them.

These commands are, for example, the loading of a specific file created by the plugin.

The final user should not see and know that the plugin calls other software and I would like that the command line of that software remains hidden.

So, how can I pass the information required by the command line directly by the python script of the plugin?

BTW, this is a code snippet within the python script that calls the external software:

if not self.chkOnlyInputMT.isChecked() and sewatFlag == 0:
        QtGui.QMessageBox.information(None, 'Information', 'Running ... \n' + mt3ddir )
        mt_namefile = transportName + '.nam_mt3dms'
        sub.Popen([mt3ddir, mt_namefile], cwd=pathfile)

In other words, I would like that, in this script, the information written in another file will be read by the command line without interacting with it..


texts: Trying to export a document 8ft x 4ft from ArcMap and it keeps coming out blank

texts:
I have a simple MXD that I have the page size properties set to 8ft x 4ft, the intended output being a large wall poster of an orthophoto.  I have set my data frame to the 8ft x 4ft size, brought in the orthophoto mosaic layer, centered it on the area of interest, and that is all, no other layers.  I then went to Export Map in the File menu and choose a destination and a name, and ArcMap then exports, but the resulting output is a blank white document.

If I change the page size to a standard large format (ie. ARCH E) and no other changes, it comes out fine with the orthophoto visible, but changing to the custom page size of 8ftx4ft I always get a blank page, though at the right size.

I have done a raster clip of the mosaic, with the clipping boundaries set as a shapefile of the 8ft x 4ft area of interest and it gave me a proper TIFF, however does someone know why I am getting a blank 8ft x 4ft page when trying to just export to PDF like a regular map?


texts:
I just tried creating the following function on the function editor on the select by expression

from qgis.core import *
from qgis.gui import *
from qgis.utils import *

@qgsfunction(args=0, group='Python')
def is_selected(values,feature,parent):
layer = qgis.utils.iface.activeLayer()
return feature.id() in layer.selectedFeaturesIds()

However, even if I press the load button, when I go back to the custom , it doesn't list it.

How to add my function?


texts:
This question may seem as it was asked before, but I think I just don't use the right terminology. Here's my problem. I'm an archaeology student trying to make an internal map of a prehistoric dwelling in which I use non-earth coordinates to show the localisation of artifacts inside the dwelling. I don't use any GPS coordinates, cause all the coordinates from my excel database are from an arbitrary x0-y0 which is a absolute point on the base drawing of the map (Is that what we call a non-earth system?).

So, I have no problem scanning the plan of the dwelling, making a png with it and using it as a raster. I put my coordinates points (8 of them) on the raster layer with the georeferencing plugin and it work. But, the problem is that I want to add vector layers of, for example, the turf ridge of the dwelling, or the rock pavement, or the different fire-pits. These vector layers are DXF that I draw on Illustrator. I'm able to send those on my qgis project, but they don't have embedded coordinate points so they appear on my base map at a totally wrong scale. My initial plan was to georeference these as you do with a raster in qgis, but I don't see how. Is there a way to resize and manually adjust those over my raster to correct the scale? My plan is to be able to work with the vector/shapefile layers and keep the scan raster invisible.


texts: Can a featureGroup be added to another featureGroup in leaflet?

texts:
I have some points that I want to be clustered in a L.MarkerClusterGroup, which works well:

var cluster = new L.MarkerClusterGroup(options);
features.foreach(function(feature) {
    var layer = L.marker(feature, {title: feature.name});
    cluster.addLayer(layer);
});

And some geofences that that I don't want clustered, in a drawable layer, which also works well:

var easel = L.FeatureGroup().addTo(map);
map.draw = new L.Control.Draw({edit: {featureGroup: easel}});

How can I add the cluster featureGroup to the editable featureGroup, so that the points get clustered and are editable?

I have tried adding the individual feature layers into both feature groups, but then I see both the cluster and all the individual points at the same time. The editing would have to occur only if you've zoomed in far enough to see an individual marker of a cluster.

Ideally, it would be great I could have a parent featureGroup, and have those two featureGroups as children then make the parent drawable.  But I haven't seen any documentation to suggest this is possible:

// psuedo code only, may not even be possible
parent = L.FeatureGroup().addTo(map);
easel = L.FeatureGroup().addTo(map);
cluster = new L.MarkerClusterGroup(options);
easel.addTo(parent);
cluster.addTo(parent);
map.draw = new L.Control.Draw({edit: {featureGroup: parent}});


texts: How to enable the snapping for a layer with the tolerance value with python programming

texts:
I have a vector polygon and vector polyline layer. I want to enable the snapping for both the layer with tolerance value? How can i do that using python programming. I have tried this code but getting the following error:

snapper = QgsMapCanvasSnapper(canvas)
snapper.setSnapSettingsForLayer(aLayer,True,SnapToVertexAndSegment,pixels,0.0001,True)

Traceback (most recent call last): File "", line 1, in  AttributeError: 'QgsMapCanvasSnapper' object has no attribute 'setSnapSettingsForLayer'

Can anyone help me in setting the snapping options to the currentLayer using python programming?


texts: How to use Google Maps API V3 to diplay circle with Openlayers

texts:
I found bellow code to display Circle by click event with Radius 60.0m.(google maps API v3). https://stackoverflow.com/questions/15111403/google-maps-api-v3-change-circle-radius-by-click-event
 I want to display a circle (within a limit of 3000m, 5000m )  and the popup window to appear attribute details of buffer zone. (school name and the school type & etc.,).
Please help me  to properly place a this JavaScript code into my  HTML page.


texts: Show QgsLayerTreeView in a custom applications

texts:
I need to read some raster data and show their information in a tree widget, but while I use QgsLayerTreeView, it can not shown in my GUI.
The last is my code, I want to know where I am wrong.

from PyQt5.QtWidgets import *
from PyQt5.QtCore import *
from qgis.gui import *
from qgis.core import *
import sys
import os

class MyWnd(QMainWindow):
    def __init__(self):
        QMainWindow.__init__(self)

        self.main_widget = QWidget()
        self.tree = QgsLayerTreeView()

        pro = QgsProject.instance()
        pro.read("D:/项目/GUI开发/code/20181030/data2/my_pro.qgs")

        k = 1
        for file in os.listdir(file_folder):
            if file[-3:] == 'tif':
                filename = file_folder + "/" + file
                file_info = QFileInfo(filename)
                base_name = file_info.baseName()

                if k < 10:
                    layer = QgsRasterLayer(filename, base_name)
                    QgsProject.instance().addMapLayer(layer)
                    k = k + 1

        root = QgsProject.instance().layerTreeRoot()
        model = QgsLayerTreeModel(root)
        self.tree.setModel(model)

        self.setCentralWidget(self.tree)
        self.show()
        pro.write()

if __name__ == '__main__':

    file_folder = "D:/项目/GUI开发/code/20181030/data2"
    QgsApplication.setPrefixPath(r"C:\OSGeo4W64", True)
    qgs = QgsApplication([], True)
    qgs.initQgis()
    app = QApplication(sys.argv)
    w = MyWnd()
    sys.exit(app.exec_())
    qgs.exitQgis()


texts: Buffer Not Taking Attribute Table Information

texts:
I have joined a spreadsheet to a polyline shapefile. The information from the spreadsheet is present in the attribute table however when I apply a buffer the information is not in the attribute table for the buffer. I have done this before and gotten it to work.

Any advice?


texts: Does GeoPackage impose a maximum number of (non-geometry) columns per feature?

texts:
Working in R, I tried to write a simple feature file with 402 rows and 1170 columns  (only one geometry column) and got the following error:

> dim(sf_object)
## [1]  402 1170

> sf::write_sf(sf_object, "sf-object.gpkg" , driver = "GPKG")

## Failed to create feature 0 in sf-object
## Failed to create feature 0 in sf-object
## Error in st_write.sf(..., quiet = quiet, delete_layer = delete_layer) :
##  failed writing to temporary file
##  C:\Users\URBAND~1\AppData\Local\Temp\RtmpukYmED\filee82cdc6ada.gpkg

When I split the simple feature object into two objects (each with half the columns) the code runs without error.

Does the GeoPackage standard impose a maximum number of columns per feature that I am unaware of?


texts: QGIS - How to globe-rotate a raster around the Greenwich Meridian?

texts:
This is probably the easiest way to explain my problem:

I'm not sure how to say what this problem is - is it that the shapefiles and the raster are using different longitudinal origins? (probably obvious I have a fairly loose grasp of GIS)

Regardless, what operation would I use to globe-rotate the raster so that it was aligned with the shape files? To be clear, it seems like the two data sets agree on coordinates and projection (WGS 84), it's just that the shape files all put the americas on the left and for whatever reason the topo raster puts them (and part of europe/africa) on the right.

I'm using QGIS 2.16.3 on a mac, and I'm fairly capable at installing & compiling code if needed.


texts: Geocoding from street layer using QGIS (MMQGIS)

texts:
I have a question about geocoding from street layer using the MMQGIS plugin in QGIS. I have a list of addresses that I would like to geocode into latitude and longitude. Because the number of records is really huge, I want to use street layer map instead of Google map or other web-based geocoding tool. However, my output did not give my what I wanted, and I would like to see if anyone can help me figure out where I did wrong.

My address file is a list of addresses with street name and house number, and it's pretty clean. I think the issue probably comes from the street line file that I used. I downloaded this shapefile from Milwaukee city government. The attribute tables looks like something like this (the street address is separated, and it has no longitude/latitude):

When using the plugin, I put "street line start" and "street line end" for X and Y attributes. Just like this:

No errors show up, but my output is empty.

Anyone knows what went wrong in this process? Maybe it's because the street name format is different in these two files? Or was it the problem with the map coordinates?


texts:
Staff has anyone ever used OpenLayers Editor?
Openlayers researching on the internet I found this link: https://github.com/geops/ole
It is an editor based on OpenLayers open source?
How do I set it on my machine to test with PostGIS and GeoServer?


texts: LayerList Legend onload visible all Layer

texts:
I am using LayerList widget and have some problem with Legend when the LayerList is loaded the first time. Actually I have a lot of Group Layers and some of them are turned off but inside the group there are some layers turned on but not visible because the Group Layer is turned off at the beginning.

The problem is the first time when LayerList is loaded, in Legend tab all symbols for all sublayers which are turned on (even on those which Group Layer is turned off and they are not visible on map) are shown.
If I turn on or off (by mouse click) some layer then the legend becomes normal and doesn't have this problem.

Did anyone have the same problem and found a solution?
I try on method onload to refresh the widget or to turn on layer and then turn off programatically but it did not solve the problem


texts:
So I've been trying to download and re-name my attachments. I have been successful on ArcMap HOWEVER, the python tool I've been able to use only downloads one attachment per datapoint. I have multiple attachments per datapoint and would like to download them all.
The code that worked for me is https://gis.stackexchange.com/a/251829/165390
How can I edit this code to download multiple attachments? Is there a better Python tool out there?


texts: Clipping Rasters from Mosaic Dataset in Original Coordinate System?

texts:
I am building a large mosaic dataset that contains several hundred rasters, which are spread across three UTM zones.  My goal is to be build a tool in ArcMap (most likely within model builder) that will clip the mosaic dataset by a user-selected shapefile, which would then be exported for use in another software program.

The secondary software utilizes 8 digit UTM coordinate systems (i.e., has the UTM zone number in front of the easting).  As such, I've made custom coordinate systems covering the study area which have the UTM zone number included in front of the false easting (e.g., UTM zone 14N would have a false easting of 14500000m and so on).

The rasters in the mosaic dataset are already projected in these custom coordinate systems, but the mosaic dataset itself is in WGS84 World Mercator, although that could be changed if it would help.

Whenever I clip the mosaic dataset via feature class, the output raster is in the coordinate system of the mosaic dataset (WGS84 World Mercator), while I would like it to be in the orginal coordinate system of the source raster.

I'm looking for a way to either retain the coordinate system of the source raster(s) when clipping or automate a tool to project them in the correct coordinate system based on their location or UTM zone.  I have added the correct UTM zone of each source raster in the mosaic dataset's attribute table, but don't know where to go from here.

I'm not very familiar with Python or programming in general, so I would prefer to create a tool in model builder, but I am not opposed to making a python script if that's the only potential method.

One idea, although I'm not sure how to execute it, would be to run an If Then type statement on the location of the feature class used for clipping.  I.e., if the feature class is within a shapefile representing Zone 15, then use 'X' projection.  If the clipping feature class is within a shapefile representing Zone 14, then use 'Y' projection, etc.

After @Michael's comments, I'll work on integrating that script.  Looking at this link http://resources.arcgis.com/en/help/main/10.1/index.html#//002w00000022000000, I am eventually trying to build a similar true/false type logic.  The goal would be to run a script testing which UTM zone the clipping shapefile is in, and then using the results of that statement to set the coordinate system for the extraction process.

Any thoughts on how to accomplish that?  I was thinking testing whether the extraction shapefile's centroid is contained within a shapefile representing an entire UTM zone.  IF I had 3 shapefiles (one for each UTM zone in the study area), the extraction shapefile could be tested somehow to see in which shapefile it is contained, thus providing a "true" response in that instance. Something along the lines of "if the extraction shapefile's centroid is within the zone 14 shapefile, use coordinate system A.  If the extraction shapefile's centroid is within the zone 15 shapefile, use coordinate system B" etc.


texts: Using simple Math in SLD file on Geoserver?

texts:
Is it possible to use simple Math in an SLD file on Geoserver? Suppose I want to read a text size property from a database and multiply its value, how can I do this?

<sld:TextSymbolizer>
    ...
    <sld:CssParameter name="font-size">
        <ogc:PropertyName>text_size</ogc:PropertyName> * 2.1
    </sld:CssParameter>
    ...
</sld:TextSymbolizer>

This simple approach obviously doesn't work, so I thought there must be some way to achieve this.


texts: Identifying high NDVI values close to vector lines

texts:
I have NDVI raster data and vector line data for the St George Queensland irrigation channels.

I am trying to detect irrigation channel leakage using high NDVI values detected on the immediate bank.

How would I identify high NDVI values close to the channel lines in ArcGIS Desktop?


texts: Adding a hotlink to an attribute field in QGIS

texts:
How to add hotlink to QGIS? I saw this but had some problem. I want to open network location like \\fileserver\filefolder\file1.excel but it ends with error.

Edit:
    The network location is stored in a column of the table eg. URL column contains URL value for each attribute. Different row has different URL.


texts:
New to leaflet. Figured a method of adding popups to show features at several geographic locations on a map. Each popup gives information via embedded URL link for a feature at that location. Displays fine. Problem arises in that ALL popups display URL of the 1st popup added to map instead of displaying their own unique URL for their individual location.

Question: Can it be rectified so that each popup displays its own assigned URL instead of adopting the URL of the first popup added to the map?

My code is as follows:

var myDataPoint = L.marker([51.491027, -0.224525]).addTo(map);

var popup = L.popup();

function onMapClick(e) {
popup
  .setLatLng(e.latlng)
  .setContent(e.latlng.toString() + ' Information about theatre"')
  .openOn(map);
}
map.on('click', onMapClick);


texts: Why doesn't QGIS 2.8.2 Road Graph accept my road network vector?

texts:
I'm using the Finnish official Department of Transport's "digiroad" road network vector file, which has all kinds of roads, cycle paths etc on it, but I can't use Road Graph on it.

If I create a new line vector and draw some lines, I can use Road Graph on it.

I don't see how the digiroad is different, other than having several different types of roads. Why can't Road Graph utilize it (it doesn't even let me choose it)?

I'm a beginner, but I've used the same material before on arcgis.


texts: QGIS - "Bad Allocation" Exception - on identify

texts:
I have QGIS (v1.8) set up to identify in Top down mode.  When I use the identify tool on the map, I frequently get a Bad Allocation exception message and I was wondering what gives rise to this error.  It seems to occur, after all layers have been scanned for intersects with the click location (i.e. the status displays as 'identifying done') - it seems to appear when the popup window is being generated.

Any thoughts on what I can try to do to resolve this?

THanks


texts:
I am writing an ArcPy script that will run from a tool with user inputs.  The user will select an input layer and output name/directory.  I have my current code below.  The tool will then do the following:

Loop through each record in the input layer and do the following within each loop cycle:

Clip a census blocks layer by the input feature
Create a field in the output clip layer to hold an area ratio value
Create a field in the output clip layer to hold an estimated population value
Perform a field calculation between two fields in the output clip layer to produce the area ratio (new area divided into original area)
Perform a field calculation bewteen two fields in the output clip layer to produce the estimated population (area ratio value calculated in #4 mulitplied by the population field)

Let's say the input layer contains three records (1500-ft, 2000-ft, and 2500-ft).  There would be three output layers that are clipped census block layers based upon the records from the input, each with estimated census populations based upon the distances from the input feature.  The layers would have names like Output_ft_1500, Output_ft_2000, and Output_ft_2500.

My code can do everything except complete the field calculations.  I'm not sure if I'm missing something simple with the syntax for a SQL expression using fields or what.

Please let me know if you have any questions. I tried to explain things as best as I could.

Here's the error I get:

Traceback (most recent call last):
  File "F:\Scripts\ArcGIS Geoprocessing\SARA Tool\EstimateCensusPopulation.py", line 73, in
    fieldExpression = newArea / orgArea
TypeError: unsupported operand type(s) for /: 'Field' and 'Field'

Failed to execute (EstimateCensusPopulation).

# Import arcpy
import arcpy

#Set workspace

# Sara Facility
sara = arcpy.GetParameterAsText(0)

# Clip Feature - U.S. Census Blocks
censusBlocks = r'\\CCPASR07\ncgs$\Scripts\ArcGIS Geoprocessing\SARA Tool\SARA_Tool_DEV.gdb\CensusBlocks_2010'

# Output
output = arcpy.GetParameterAsText(1)

# Search cursor
cursor = arcpy.SearchCursor(sara)

for row in cursor:
    # Clip features
    feat = row.Shape

    # Append buffer distance and units to name
    # Buffer distance
    buffDist = str(int(row.BUFFDIST))

    # Buffer units to name
    buffUnits = row.UNITS

    # Appended output name variable
    buffAppend = '_' + buffUnits + '_' + buffDist

    # Execute clip tool on each row
    newInput = arcpy.Clip_analysis(censusBlocks, feat, output + buffAppend)

    # Add message that Clip is completed
    arcpy.AddMessage('Feature Clip operation completed')

    # Add field to hold clip area to original area ratio
    areaRatioFieldName = 'AREARATIO'
    areaRatioFieldType = 'DOUBLE'
    # Execut tool
    arcpy.AddField_management(newInput, areaRatioFieldName, areaRatioFieldType)

    # Add message that Area Ratio Field has been added
    arcpy.AddMessage('Area Ratio field added')

    # Add field to hold estimated population
    estPopFieldName = 'ESTPOP'
    estPopFieldType = 'LONG'
    # Execut tool
    arcpy.AddField_management(newInput, estPopFieldName, estPopFieldType)

    # Add message that Estimated Population Field has been added
    arcpy.AddMessage('Estimated Population field added')

    ### Code does not work after this point ###
    areaInField = arcpy.ListFields(newInput, 'AREARATIO')[0]
    newArea = arcpy.ListFields(newInput, 'Shape_Area')[0]
    orgArea = arcpy.ListFields(newInput, 'ORAREA')[0]

    fieldExpression = newArea / orgArea
    arcpy.CalculateField_management(newInput, areaInField, fieldExpression, 'PYTHON_9.3')

del cursor


texts: How to change value of an individual pixel using arcpy?

texts:
I want to change the value of an individual pixel in a raster using arcpy. The code I wrote for changing the value is below, but it changes all of pixels in my raster. How can I set an individual pixel in my code?

import arcpy
import numpy
raster = arcpy.Raster(r"C:\test.jpg")
array = arcpy.RasterToNumPyArray(raster)
(height, width)=array.shape
for row in range(0,1):
    for col in range(0,1):

        newarray = array * 2.0

        newraster  = arcpy.NumPyArrayToRaster(newarray)
        newraster.save(r"C:\export.gdb\t")

Thanks


texts:
I try to open the attribute table of a layer witch has a PythonInit fonction with a Ui form and returns "AttributeError: 'NoneType' object has no attribute 'length'".
Made so far:

from PyQt4.QtCore import *
from PyQt4.QtGui import *

amonte = None
aval = None
panta = None
layer = None
lungimeTextEdit = None
ValLungime = None
valoare = None
feature = None

def formOpen(dialog,layer,feature):
   global myDialog
   myDialog = dialog
   global lungimeTextEdit
   global amonte
   global aval
   global panta

   lungimeTextEdit = dialog.findChild(QLineEdit,"shape_leng")
   amonte = dialog.findChild(QLineEdit, "amonte")
   aval = dialog.findChild(QLineEdit,"aval" )
   panta = dialog.findChild(QLineEdit, "pente")
   buttonBox = dialog.findChild(QDialogButtonBox,"buttonBox")
   lungimeTextEdit.setText(str(feature.geometry().length()))
   buttonBox.accepted.disconnect(myDialog.accept)


texts: meeting point minimizing travel distance for participants?

texts:
The problem is to find the point minimizing the travel distance for around 100 persons in different regions who want to meet in the same place. Travel is by car not by plane.

Assuming that I get access to an API giving me mileage / kilometric distance in terms of highway travel between any two points, how can I find the best place to meet?


texts: Loading color table (.csv) file to interpret raster data using C#?

texts:
I have a raster (GeoTIFF .tif file) that is stored in three bands (Red, Green, Blue), with each band's data type as a byte. I have a separate color table in the form of a .csv file, that stores a height value (in meters), a Red value, Green value, and Blue value on each line. I need to able to utilize this color table to interpret the RGB values in the raster as raw height values. I am using C# bindings of GDAL, OGR and OSR. How do I get the file loaded so that it can be represented as a GDAL color table?

I know OGR has a CSV driver which you can use to Open the file, which I have done. This creates a DataSource object. But I have no idea how to convert this DataSource to a GDAL ColorTable (if it is even possible). Or is another method needed, such as reading the file line by line manually?


texts: How to change the shape of a multidimensional raster?

texts:
I have a hyper-spectral data file which is acquired with pushbroom method. This means basically data cube has the shape like (swath, spectral bands, along the track). What I mean is the shape differs from the common shape(spatial x, spatial y, spectral bands) format.
I can use this data on python with numpy, simply by exchanging the axes. But I need to perform some experiment on QGIS, which I don't have any familiarity.
Is it possible to perform this on QGIS, and if yes, how to do that?


texts: How to Mosaic DTED files (RasterDataset) using C#?

texts:
I use ArcGIS .NET SDK in C#.

I'm able to load some DTED files in my ArcMap and show them in proper colors. What I want is to concat/merge/mosaic some DTED files. How can I do this?

I want to mosaic IRaster or IRasterDataset. I googled A LOT and struggled with ArcObjects Help, but still unsuccessful.


texts: This Program/Project is working. Using the centroid function, I am able to filter the layer but I am confused about how to add another function where I can get the integer value from the column "active_cases".


texts: GeoServer 2.10.4 - Vector layer always on top

texts:
I have successfully created a new vector layer in my GeoServer 2.10.4. It shall be displayed on top of a basemap and an additional raster layer, but it displays below the raster. I want this order:

Basemap
Raster
Vector

But got this one:

Basemap
Vector
Raster

How can I change this, so the vector layer always sits on top of everything else? From various similar questions I found that it should be achievable using the SLD of the layer and inserting z-order=X somewhere, so here is my SLD for this particular layer. I tried values of 0, 1 and 99 so far, all with the same result.

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor version="1.0.0"
 xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
 xmlns:ogc="http://www.opengis.net/ogc"
 xmlns:xlink="http://www.w3.org/1999/xlink"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <!-- a Named Layer is the basic building block of an SLD document -->
  <NamedLayer>
    <Name>default_polygon</Name>
    <UserStyle>
    <!-- Styles can have names, titles and abstracts -->
      <Title>Default Polygon</Title>
      <Abstract>A sample style that draws a polygon</Abstract>
      <!-- FeatureTypeStyles describe how to render different features -->
      <!-- A FeatureTypeStyle for rendering polygons -->
      <FeatureTypeStyle>
        <Rule>
          <Name>rule1</Name>
          <Title>Red Polygon with Black Outline</Title>
          <Abstract>A polygon with a red fill and a 1 pixel black outline</Abstract>
          <PolygonSymbolizer>
            <Fill>
              <CssParameter name="fill">#ff0000</CssParameter>
            </Fill>
            <Stroke>
              <CssParameter name="stroke">#000000</CssParameter>
              <CssParameter name="stroke-width">1</CssParameter>
              <CssParameter name="z-index">99</CssParameter>
            </Stroke>
          </PolygonSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts:
I know that I can just call one of any number of geocoding APIs to perform this for me. This is more of a curiosity than anything else.

With that aside, how does someone (like Google) turn a city name into a latitude and longitude? From playing around with Google Maps it appears that most of the time a city is pinpointed to it's downtown area, but there seems to be no reason why one street corner is picked over any other street corner.

I am just curious what algorithm a geocoding service provider would utilize for performing such an action?


texts:
I am doing an ogr2ogr process where i will be making a join between two tables, imported as shapefiles and .csv to join on a column with idential attributes.

The code i've made so far looks like this:

    ogr2ogr -sql SELECT spmgnt100all_netto_rumnr.*, Areas_editeret.csv FROM spmgntall_netto_rumnr.* LEFT JOIN Areas_editeret.csv ON spmgnt100all_netto_rumnr.RUMNR = spmgnt100all_netto_rumnr.RUMNR Areas_Joined.shp spgmgnt_Joined.shp

But I am getting an error: "Error 1: Failed to create directory spgmnt100all_netto_rumnr., for shapefile datastore"
ESRI Shapefile driver failed to create spgmnt100all_netto_rumnr.,

I suppose the error occurs when the .shp-table is submitted into another format and that is where the driver fails?
Or maybe it means that ogr2ogr tries to make a shapefile that is the result of the process. And which fails.

But i might also understand the syntax wrong. I've tried to read my way through http://www.gdal.org/ogr_sql.html


texts: Drawing arc between two points using Java?

texts:
I have two points (x1y1 and x2y2), the centre point between the two lines (x0y0), and I am trying to draw an arc using the following Java method:

g.fillArc(x of initial point, y of initial point, width of arc, height of arc, start angle, end angle)

to get this:

I have found radius:

int r = (int)Math.sqrt((x1-x0)*(x1-x0) + (y1-y0)*(y1-y0));

then width and height:

int width = 2*r;
int height = 2*r;

and then start and end angles:

int startAngle1 = (int) (180/Math.PI*Math.atan2(y1-y0, x1-x0));
int endAngle1 = (int) (180/Math.PI*Math.atan2(y2-y0, x2-x0));

Instead, I get this:

or this:

Trying to figure out the proper way to draw the arc.

EDIT:

These are geographical coordinates as well as x and y coordinates on the screen based on the zoom property:

Point 1:

X1 coordinates: [25.7745, -81.4589]
X1 and Y1: 574043 : 893092
X2 coordinates: [25.82111, -81.45741]
X2 and Y2: 574052 : 892790
X0 coordinates: [25.797805001897803, -81.45815514647514]
X0 and Y0: 574047 : 892941

Point 2:

X1 coordinates: [26.85249, -80.8777]
X1 and Y1: 577429 : 886086
X2 coordinates: [26.87219, -80.8392]
X2 and Y2: 577653 : 885957
X0 coordinates: [26.862341303504717, -80.8584516762013]
X0 and Y0: 577541 : 886022


texts: How many of you out there work from home? What kind of company / organization do you work for (private, government, municipal, institution..?). How did you convince your employer it was a good idea? What are the major drawbacks to working in a GIS position at home?


texts:
As far as I know Geoserver web administration page language depends on the language settings in browser. So if my preferred language is Russian I see partially translated page, and if I need it in English, I just set English before Russian in browser settings.
Is there any way to override this behavior by editing some xml config or properties file? I need admin page to be always displayed in English, regardless of user browser language.
I am using Geoserver 2.6.1.


texts: Writing selected features to shapefile using OGR and Python?

texts:
I have two shapefiles that intersect: one containing a line, and the other containing polygons. I'm trying to use OGR and Python to create a buffer around the line and select polygons within it, then write those features to a shapefile.

##import modules
import ogr
import gdalconst
import sys
import os

##Define input files
power = (r'C:\Users\Megan\Desktop\powerlines\PowerLine.shp')
parcel = (r'C:\Users\Megan\Desktop\parcels\Parcels.shp')
outName = ("BufferedParcels.shp")

##Define buffer width
bufferWidth = 250

##Assign driver
driver = ogr.GetDriverByName('ESRI Shapefile')

##Open files
powerData = driver.Open(power, gdalconst.GA_ReadOnly)
parcelData = driver.Open(parcel, gdalconst.GA_ReadOnly)

##Verify files were opened correctly, exit if not
if powerData is None:
    print('Could not open ' + power)
    sys.exit(1)
if parcelData is None:
    print('Could not open ' + parcel)
    sys.exit(1)

inputDS = driver.Open(power)
if inputDS is None:
    print('Could not open', power)
    sys.exit(1)

##Get feature layer and geometry from powerlines
powerLayer = powerData.GetLayer(0)
powerFeature = powerLayer.GetNextFeature()
powerGeom = powerFeature.GetGeometryRef()

##Get feature layer and count from parcels
parcelLayer = parcelData.GetLayer(0)
parcelFeatureCo = parcelLayer.GetFeatureCount()

##Re-write output file if it already exists
inputLayer = inputDS.GetLayer()
if os.path.exists(outName):
    os.remove(outName)

##Apply buffer to powerline geometry, give affordance that it worked
powerBuffer = powerGeom.Buffer(bufferWidth)
print('Created 250 ft buffer around',power)
##Attempt to get output driver, exit in case of failure
try:
    outputDS = driver.CreateDataSource(outName)
except:
    print('Could not create',outName)
    sys.exit(1)

##Define new layer, exit in case of failure
newLayer = outputDS.CreateLayer('BufferedParcels',geom_type = ogr.wkbPolygon)
if newLayer is None:
    print('Could not create layer for buffer in output DS')
    sys.exit(1)
##Get new layer definition
newLayerDef = newLayer.GetLayerDefn()

##Give affordance that program is running
##Set feature ID
##Loop through features in parcelLayer
print ('Writing the following parcels to shapefile:')
featureID = 0
for i in parcelLayer:
    ##parcelFeatures = parcelLayer.GetFeature(i)
    parcelGeom = i.GetGeometryRef()
    bufferContains = powerBuffer.Contains(parcelGeom)

    names = i.GetField('SITUSADDR')
    area = i.GetField('AREA')
    relevantVals = (names, area)
    fmt = '%s %s'
    ##Determine if features are within buffer
    ##Print stats if 'True'
    if bufferContains == True:
        print (fmt % relevantVals)
        try:
            newFeature = ogr.Feature(newLayerDef)
            newFeature.SetGeometry(poly)
            newFeature.SetFID(featureID)
            newLayer.CreateFeature(newFeature)
        except:
            ##print('Error printing shapefile.')
            newFeature.Destroy()
            ##parcelLayer.Destroy()
        featureID += 1
    i.Destroy()

##Give affordance that program worked
print('Wrote', featureID, 'parcels to shapefile.')

inputDS.Destroy()
outputDS.Destroy()

The code above successfully writes a list of the polygons I'm looking for and creates a shapefile. However, the file has an empty attribute table and doesn't show any geometry.

Can anyone see where I'm going wrong?


texts: Can Trimble Business Center export projection information in any useable GIS format?

texts:
In Trimble Business Center, a surveyor can create a local ground system for a given location easily. For example, when building a bridge, the surveyor can create a ground system based on State Plane so that an engineer can be sure that their design can be built correctly. Later in the project, I'm asked to take external data and project it to this local coordinate system.
Can TBC export a .prj file? If not, can it export any other representation of this custom ground such as WKT that I can then use to create a .prj file, and then project any data into the projects CS?
Currently, my surveyor is only able to show me a screenshot showing the False Easting, False Northing, and scale factor for the local CS. I can incorporate these into a prj file myself, but it is a bit time consuming and prone to errors. I'd like a reliable way for a surveyor to create a custom projection that I can easily use in a GIS.


texts: How to print Large map in QGIS in normal paper size(A4,LEGAL)?

texts: How to print large map in A4 size like "Tile printing"method where you can automatically divide large map into equal parts(A4) then join it later by pasting, i tried ATLAS plugin but it gives me different output, it divides the large map into A4 size but it has big overlaps in every A4 paper it produced.


texts: I am new to Mapbox studio so I may be overlooking something very basic - I have a polygon shapefile that I am trying to style - I want one label per polygon but am getting more - often four labels.  I have tried using the buffer-size and text-repeat-distance options but I can't seem to get it to work.


texts: How do you retrieve an accidentally deleted point in an ArcGIS Web Mapping Application?

texts:
Hypothetical scenario: I've shared an ArcGIS Web Mapping Application with non-tech-savvy staff members. One of the staff members accidentally deletes a point within the Web Mapping Application.

Is it possible to undo this operation within the Web Mapping Application?
If not, Can I retrieve the data from the linked database?


texts: Problem to transform events from one route to another

texts:
I have a road feature class layer containing parts like 1_10, 1_20. 2_10, 2_20... etc. 1_, 2_ is later used to create the route and _10, _20 is indicating the direction (the road feature classes were digitized in the same direction - they told me at least!).
In this feature classe I calculate the "from" and "to" values for each segment.
Ater that I create a route with arcpy.CreateRoutes_lr that seems to work great.
Now for this route, linear events are collected over time.

After some time there are changes in the road feature class.
So my arcpy script now recalculates the new lengths of the road segments and their "from" and  "to" values, then creates a new route. This all seems to work fine BUT I am not able to transform the old events from the old route to the new route with arcpy.TransformRouteEvents_lr. I get a lot of event artifacts (like 0.5m long) and even events from lets say route 7 are "mapped" to route 10, or Event X is on route 10 and 11... However I think the most errors are like ths (the events are basing on the original segments):

Event1 on route1 from 0 to 20
Event2 on route1 from 20 to 30
Event3 on route1 from 0 to 25 ??!!!

Another problem is that my routes consists of loops which I think is an additional challenge.

I also tried a workaround as creating a new feature class with the old routes and the old events and then with arcpy.LocateFeaturesAlongRoutes_lr to map it on the new routes but that does not seem to work either (actually I have quite similar results as with the transform function).

I dont want to map old events on a changed/new routes (since that makes no sence for pavement). And this seems to work at least. The problems are events on routes that have not changed, but maybe the from/to values of the events must be recalculated.

I hope this wall of text makes sense somehow and unfortunately I can not post pictures yet.
Any help is really appreciated!

Edit: after doing some research I found out that all features that are not a single line create errors so for example Y-or T-junctions in a route, or looping routes.


texts: Building Graph in QGIS for Service Area computing

texts:
I do have

a vector shapefile with a streets network
a point shapefile with 5 supermarkets as points

Everything is right now in WGS 84.

I would like to perform QNEAT's "ISO AREA AS POLYGON from Layer" around my Supermarkets.

First, I think I have to build a graph from my streets Layer. But that requires a projected CRS such as UTM 32 N?
How to transform WGS correctly in UTM 32 N under QGIS? (Got only failures)

Then I have to perform "Build Graph" on the Streets UTM32N? Right?
And then I can perform ISO Area on both shapefiles.
The supermarkets must be in UTM 32 N as well, but otherwise, they can stay as a vector layer, is not it?


texts: Viewing ArcGIS Pro layer as JSON to illustrate CIM (Cartographic Information Model)

texts:
In order to work with the Cartographic Information Model (CIM) from ArcPy with ArcGIS Pro I am interested in viewing the JSON for a particular layer.

In a YouTube video from Esri's 2019 Developer Summit:

  Product Engineer, Jeff Barrette, demonstrates new and improved tools
  for map creation automation, to include features: advanced symbology,
  arcpy.mp, and the Cartographic Information Model (CIM).

At the 2:58 mark of that video, the presenter shows a JSON file for an ArcGIS Pro layer as background for working with the CIM.

How was such a JSON file obtained?


texts:
I have a fruit orchard dataset. One of the columns contain the year the orchard was planted. In  another column I would like to have the orchard age. However this pose a problem, since the age changes every year.

My question:
Is it possible to auto calculate orchard age based on the year column in postgis? If it is possible, please specify how.


texts: Format multiple values in a label with different colors in QGIS

texts:
I need to create complex labels for sites that contain multiple values in each. These values are colored depending on their value (ex. in range = green, out of range = red). I need to do this from an expression. Has QGIS implemented this feature yet?

Here's and image of what I'm trying to reproduce:


texts:
I have single band DEMs I'm trying to use in an ImageMosaic in GeoServer 2.11.2 that displays a black background color by default. I want that background color to go away and the default color ramping to range of displayed values to occur.

The default behavior is almost exactly what I want - when you zoom in, it scales the colors to the min/max of the dataset currently in view.

All I want to do is hide the background color of -9999.0 so that scaling isn't messed up by the nodata color. However all attempts to fix this don't seem to work.

-Setting InputTransparentColor and OutputTransparentColor cause rendering to fail with an error 'Unsupported data type'.

-Setting an SLD ColorMap that hides -9999.0 values 'works' but you need to then to specify the rest of the color ramp which I don't want to do because it needs to be dynamic. (So I don't want this answer).

-Setting an SLD ContrastEnhancement causes a NullPointerException

<ContrastEnhancement>
  <Normalize />
</ContrastEnhancement>

And trying to set VendorOptions inside the ContrastEnhancement Node doesn't validate.


texts: Adding ArcGIS Server feature service to ArcMap data frame using Python?

texts:
There's an external ArcGIS Server that I've been manually going into through ArcCatalog to pull down feature services and save them off as feature classes. I do this twice a month, so I would like to explore if it's possible to do this automatically with Python/ArcPy.

Is this possible?

Current process is that I save off 4 feature classes from 3 feature services. My process is to make a connection to the server in Catalog --> drag the FS from Catalog into my TOC --> right click on the layer of interest and save off as a feature class to my local GDB.


texts:
Is there a way to loop through every selected layer in qgis and save it as a .tif file with the same layer name using a python script?

This is the code I am using, but it does nothing at all. It runs but does nothing.

Can someone help me?

from PyQt4.QtCore import *
from qgis.core import *
from qgis.gui import *
import os
import sys

layers = iface.layerTreeView().selectedLayers()
for layer in layers:
    output='Desktop'%layer.name%'.tif'
shareeditundeleteflag


texts: Display a map using geotools java in a javafx application

texts:
i have added geotools java to my Netbeans IDE and now i want to display this map in my javafx application.I have downloaded its geotiff from the same site and it's a zip file (176 mb)  which contains three files:

G8241_E1_1846_C3.j2w    ( 1 kb)
G8241_E1_1846_C3.jp2    ( 175788  kb)
G8241_E1_1846_C3.xml    (13  kb)

I have no idea how to approach that.


texts:
I am storing OS Open Local data, as downloaded here: https://www.ordnancesurvey.co.uk/opendatadownload/products.html (top result).

I've imported the .shp files.

Running the .prj through here: http://prj2epsg.org/search produces:

    PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["D_OSGB_1936",
SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",-100000.0],
PARAMETER["Central_Meridian",-2.0],PARAMETER["Scale_Factor",0.9996012717],
PARAMETER["Latitude_Of_Origin",49.0],UNIT["Meter",1.0],AUTHORITY["EPSG",27700]]

So in my postgis db I run:

SELECT srid FROM spatial_ref_sys WHERE srtext LIKE '%British_National_Grid%';

Which produces two results:
27700 and 7405

By inspection the former corresponds exactly to my .prj

Verifiying my SRID on my column:

SELECT Find_SRID('', 'buildings', 'geometry');

 find_srid
-----------
         0
(1 row)

Fine, so casting the column to 27700 should do the trick:

SELECT *
FROM buildings AS b
WHERE ST_DWithin(ST_SetSRID(b.geometry, 27700),
ST_Transform(ST_SetSRID(ST_MakePoint(-0.1047, 51.5236), 4326),27700),
500);

(Note that (51.5246, -0.1047) is downtown London)

This query produces no results. interestingly, increasing the radius to 15000, produces 8748 results. Casting an arbitrary one to 4326 and using ST_AsGeoJSON to produce a readable result gives a location in the suburbs (northwest) of London which suggests that my coordinates are somewhat off...

EDIT: Flipping the coordinates to correspond to more traditional coordinate systems also generates no results (more or less at any radius, I've gone up to 500000):

SELECT *
FROM buildings AS b
WHERE ST_DWithin(ST_SetSRID(b.geometry, 27700),
ST_Transform(ST_SetSRID(ST_MakePoint(51.5236, -0.1047), 4326),27700),
500000);

How can I figure out where and/or how things are going wrong? And how can then fix this?


texts: Is it possible to have tool parameters appear only if clause met?

texts:
I'm writing an arcpy script for a tool that has a large number of parameters (around 10).

About half of these are "dependent" (perhaps not in the technical sense of the term), as in [first parameter] "1. Add a new layer or use existing" (binary 1 for new layer or 0 for old layer possible response)
[next parameter] "2. {if 1 in first parameter} Specify workspace for new layer

It is in this way that I am describing the parameters as dependent.

My question is whether it is possible to have parameter 2 appear only if parameter 1 has a response that warrants it? I.e. if the user uses the layer in the existing .mxd then there is no need to specify the workspace.

I've looked through the various parameter options, but there doesn't seem to be a way to do this. Any advice would be greatly appreciated


texts: Same ID to selected features in different layers in QGIS

texts:
I have two different datasets that represent buildings: OSM (polygon) and a reference dataset (both point and polygon). I need to match the same buildings in these datasets, and see how these two data sets match: (number of 1-1, 1-0, 0-1, 1-n, m-1, m-n matches for polygons and points). In order to achieve this, I need to provide an ID to all of the matching features.
An example scenario is shown below (OSM represented by orange polgons and reference dataset in purple):

Is it possible to select matching buildings in different layers (manually) and automatically label them?
For example, to assign IDs in situation A, I will select all of the three features, and these three features will have the same unique ID. It is even difficult to select features belonging to different layers at the same time, since only the topmost layer's features are selected.


texts: Run QGIS script without going through QGIS

texts:
I have the following script that I'd like to run without using QGGIS. Is there a way to do that?
import processing
from processing.core.Processing import Processing
Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

input_point=r'C:\Users\Dylan\Desktop\Testing\Points2lines\Tracks.shp'
output_line=r'C:\Users\Dylan\Desktop\Testing\Points2lines\TracksLines.shp'
processing.run("qgis:pointstopath",{'INPUT':input_point,'ORDER_FIELD':'DateTime','OUTPUT':output_line, 'GROUP_FIELD':'RUNLINE_ID'})

When I try to run it I have the following error: ModuleNotFoundError: No module named 'processing'
Thnaks in advance


texts: Joining (concat) list of similar dataframes in geopandas?

texts:
I have a list of dataframes, all of the same structure and all are polygons.

How can I add them all together using geopandas?


texts: Error when Creating Layer Package in Arcpy

texts:
I'm trying to create a layer package from a feature layer (whose data source is a feature class) using arcpy. However, I keep getting the following error:

Object: Error in executing tool

Here's the code I'm using:

import arcpy

arcpy.env.overwriteOutput = True
arcpy.env.workspace = "D:/Data/"
mxd = arcpy.mapping.MapDocument("CURRENT")
print "environments set"

lpk = "layerPackage.lpk"
summary = "A summary"
tags = "tag1", "tag2"
print "layer package variables defined"

dataFrame = arcpy.mapping.ListDataFrames(mxd, "*")[0]
layers = arcpy.mapping.ListLayers(mxd, "", dataFrame)
print layers

for lyr in layers:
    try:
        if lyr.name == "Layer_I_Want":
            print "converting Layer_I_Want to lpk"
            arcpy.PackageLayer_management(lyr, lpk, "CONVERT", "PRESERVE_ARCSDE",
            "#", "ALL", "ALL", "CURRENT", "#", summary, tags)
            print "lpk created"
    except Exception as e: print(e)

mxd.save()

How do I fix it?

I've tried different combinations of "CONVERT"/"PRESERVE" but still get the same error.


texts: I have a vector map that has two lines. I want to have the two connected. Is there a way for me to draw a line to connect the two? I am using GRASS gis by the way. The two lines that I want to connect are shown here:


texts: How to extract all lines that do not intersect with buffer?

texts:
I am working in ArcGIS Desktop 10.0 and/or 10.1.

I have two road network feature classes in a file geodatabase. I have created a 5 meter buffer around the first one, so it is now a polygon feature. The other feature is a line feature. I want to intersect the line feature road network with the buffered one and generate an output feature class that includes only the road lines that fall outside of the buffer.

When I run Intersect, the output is only the portion of the lines that fall inside the buffer.

Is there a process to "clip off" the lines that fall inside the buffer polygons and just output a line feature containing all the lines outside the buffer?

I am open to using other software if necessary. Thanks for your help!


texts: I have two Layers which I added to the map using WMS service. Is there a way to perform spatial queries using Lealflet, for example; Intersection, Distance...?


texts: Can someone explain how to do advanced sorting of an attribute table using QGIS?. Every time I try to sort by script, it does not give me what I need and I don't know if there is any plugin that can do that.


texts: Extract DEM and DTM information from LIDAR image using E-Cognition

texts: I´m trying to extract DEM and DTM from a LIDAR file. Currently, I rasterize the image but don´t know which algorithm I need to apply to obtain the DEM and DTM information in different models to be able to compare both of them.


texts: Z-index in OpenLayers 3: layer ordering in OL3

texts:
Is there a method to change Z-indexes of layers in OpenLayers3 like it was in old version?

map.setLayerIndex(markers, 99); //set the marker layer to an arbitrarily high layer index

I need to change order of layers throughout using a map. So such possibility as defining z-index like this do not help

var geoLayer = new ol.layer.Vector({
    source : new ol.source.GeoJSON({
        projection : 'EPSG:900913',
        url : './myGeoJson.json'
    }),
    style : function(feature, resolution) {
        var text = resolution < 5000 ? feature.get('name') : '';
        if (!styleCache[text]) {
            styleCache[text] = [new ol.style.Style({
                fill : new ol.style.Fill({
                    color : 'rgba(255, 255, 255, 0.1)'
                }),
                stroke : new ol.style.Stroke({
                    color : '#319FD3',
                    width : 1
                }),
                text : new ol.style.Text({
                    font : '12px Calibri,sans-serif',
                    text : text,
                    fill : new ol.style.Fill({
                        color : '#000'
                    }),
                    stroke : new ol.style.Stroke({
                        color : '#fff',
                        width : 3
                    })
                }),
                zIndex : 999
            })];
        }

    }
});


texts:
I am able to set data source encoding through this answer to Setting character encoding for all layers at once using PyQGIS? and now i want to list all the datasource encodig types. how it is possible i have also check official API but nothing found!

for layer in QgsMapLayerRegistry.instance().mapLayers().values():
    layer.setProviderEncoding(u'UTF-8')
    layer.dataProvider().setEncoding(u'UTF-8')
    print layer.name(), layer.dataProvider().encoding()


texts: Prevent PostGIS Topology from splitting edges on line insert

texts: I'm dealing with a large amount of road data represented with PostGIS Topology types. Anytime I add an edge over another edge, the edge is automatically split and a new node is created at the intersection of the two lines. This creates a problem when I have roads that overlap other roads at bridges or underpasses. Is it possible to prevent edges from automatically being split when inserting an edge using TopoGeo_AddLineString? Or is this an incorrect usage of Topology?


texts: Export to TAB with styling using ogr2ogr

texts:
We have a lot of our data stored in a PostGIS database and I need to create some styled MapInfo TAB file. So I have added an MI_STYLE column (varchar254) that contains a simple brush fill

BRUSH(fc:#0000FF);PEN(c:#000000)

I then export to TAB using ogr2ogr and open in MapInfo and I see polygons with the default black line but it should be a blue polygon with a black line.

Having read this similar question it seems I should be using OGR_STYLE as I am doing styling per row. So I changed my PostGIS database column name and did the export again. Unfortunately I still get the same result as above.

Can anyone help with an example?

What column name is MapInfo expecting?

And what is the correct ogr2ogr export?


texts:
Using the Cesium library, I want to find shadows that objects fall onto each other. With Cesium, I can create shadows of objects, but I have no idea how to get the shadows that fall on objects.

How can I obtain the areas of the resulting shadows? How can I calculate the amount of shadows or amount of light falling on objects?


texts: Normalize XLabel on NDVI Index chart in Google Earth Engine

texts:
I'm using the "ui.Chart.Image.Regions" type to produce three graphics for three different SENTINEL-2 images, with EVI index values.

I would like to normalize the horizontal label (the X label) so they look the same on all three graphs.

I used the "xLabels" property of this chart type by putting the values xLabels = [0.0, 0.5, 1.0],

however I got the following error message:

  Error generating chart: The length of the xLabels list should equal
  the number of bands in the image. Found 1 bands and 3 elements in
  xLabels. ".

Here is my code: EVI INDEX


texts: Determining if line is snapped to nearest linestring in PyQGIS?

texts:
I have a polyline layer and I want to find if a linestring is snapped to the nearest linestring or not in QGIS. How can I do that in QGIS?

I am checking for the intersection but the intersection code is not working because the line string is snapped. So I want to find the snapping linestrings.

def lines(l_layer):
    for ln1,ln2 in itertools.combinations(l_layer.getFeatures(),r=2):
        print ln1.id(),ln2.id()
        if ln1.geometry().intersects(ln2.geometry()):
            print('intersects')
            print ln1.id(),ln2.id()
        else:
            print('not intersects')

Here I am not getting one linestring which is touching with the other because it is snapped and it doesn't have the vertex. So I want to find if the line is snapped.


texts: ModelBuilder routine to remove spatially duplicate features

texts:
I'm attempting to build a ModelBuilder routine in ArcGIS 10.1 based off of an existing post on Stack Exchange.  The response from 'Get Spatial' gave a great screen shot of a model that would run an input through an interative process to remove spatially duplicate features.  I've attempted to build this model on my own but have run into two problems.

I have been unable to get the 'Output Feature Class' at the end of the routine to recognize the Input Features as a Feedback Variable.  I've attempted making this connection through the Class Properties of the Output and the connect tool by picking the output and then the input.

What sort of restrictions on the Feedback Variable could be preventing me from making this connection?

I've attempted to run the model several times without the Feedback Variable connected and gotten an error that there is no input or fields that are required.  I've entered this information in the parameters.

What reasons would cause this error to occur?  Could it be related to the missing feedback Variable?

Appreciate the help.  I may need to reconsider the use of the Feedback Variable and the design of the model if I can't get around these issues.  Thanks,
Brian


texts: Is there a public OGC Catalog Service (CSW) available?

texts:
Is there an "almost always available" service URL for a public Catalog Service (CSW) available anywhere?

I'm keen to test CSW client capabilities without having to configure a CSW server.


texts: ggplot gridded SpatialPolygonsDataFrame producing odd triangles and subsetting data based on point data

texts:
I have cross posted this on SO.
Using the code below I can plot the following:
This code is adapted from here

As you can see there are  few issues with the plot. I am struggling to

Remove weird lines in plot
Only plot cells (grids) where there are data
Plot ID (see gridSpatialPolygons$values) on top of the grid cell

I realise there are a few points to this question but I hope one solution solves all.
# Load libraries
library(sp)
library(raster)
library(ggplot2)

# Projection
wgs.84 <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

# Load data
x <- c(76.82973, 76.82972, 76.82969, 76.83076, 76.83075, 76.83071, 76.83129, 76.83126, 76.83125)
y <- c(28.26734, 28.26644, 28.26508, 28.26778, 28.26733, 28.26507, 28.26912, 28.26732, 28.26687)
z <- c(-56.7879, -58.22462, -58.4211, -55.75333, -58.55153, -56.38619, -56.11011, -58.17415, -59.77212)

# Create data frame
dataset <- data.frame("LONGITUDE" = x, "LATITUDE" = y, "VALUES" = z)

# Create SpatialPointsDataFrame object
datasetSP <- SpatialPointsDataFrame(coords = dataset[,c(1,2)], data = data.frame("id" = 1:nrow(dataset), "values" = dataset$VALUES), proj4string = wgs.84)

# Extent
extentDatasetSP <-extent(datasetSP)

# Make grid options

# Cell size (map units)
xCellSizeGrid <- 0.001
yCellSizeGrid <- 0.001

# Grid
grid <- GridTopology(cellcentre.offset = c(extentDatasetSP@xmin, extentDatasetSP@ymin),
                     cellsize = c(xCellSizeGrid, yCellSizeGrid),
                     cells.dim = c(3, 7))

# Create SpatialGrid object
gridSpatial <- SpatialGrid(grid = grid, proj4string = wgs.84)

# Convert to SpatialPixels object
gridSpatialPixels <- as(gridSpatial, "SpatialPixels")

# Convert to SpatialPolygons object
gridSpatialPolygons <- as(gridSpatialPixels, "SpatialPolygons")

# Add 'id' and 'values' to every polygon
gridSpatialPolygons$id <- 1:nrow(coordinates(gridSpatialPolygons))
gridSpatialPolygons$values <- paste("Gridvalue", 1:nrow(coordinates(gridSpatialPolygons)), sep = ":")

# Get attributes from polygons
samplePointsInPolygons2 <- datasetSP %over% gridSpatialPolygons

ggplot(gridSpatialPolygons, aes(x = long, y = lat)) +
  geom_polygon(color = "red") +
  geom_point(data = dataset,
             aes(x = LONGITUDE,
                 y = LATITUDE))


texts: Converting Shapefiles to KML and extruding on height

texts:
I am working on a free/open source workflow to ingest OSM data and convert it to a KML so I can view it in 3D with Google Earth.

I used ogr2ogr to convert to a shapefile and have extracted the height values into their own field in the file. I now want to use ogr2ogr to convert to KML. I understand that I need to use the LIBKML library and extrude on my height field, but while I try to do this I repeatedly get flat KML output. I tried to look through the GDAL documentation, but there is little info on how to extrude height from a field so it shows up in the resulting KML.

This is the command I attempted:

ogr2ogr -f libkml "E:/filepath.kml" "E:/filepath.shp" -zfield "height" -t_srs "EPSG:4326"

I would really like to use ogr2ogr to solve this, but if there's another open source library I'm open to that too.


texts:
I am working with a set of coordinates from a csv file that creates a grid (of points) in QGIS. There is an id attribute along with the X and Y coordinates for each point.

I need to work with a polygon grid made out of these points. Using the grid tool doesn't work as there is a slight rotation to the points:

I need to create a grid that has a point in each intersection of lines and would look like the following grid I drew by hand:

I looked everywhere for a way of doing so, I found this tutorial but it didn't work for me.


texts: Separating two label names in QGIS in order of having a better display

texts:
I am working with layers that have several characteristics in the attributable, and I need to show two of them together.
Here is a demo map of the US state that shows me problem.

I used the string connection functions in the expression dialog for the simultaneous displaying of the "postal" and the "name" labels.

Now this two labels are displaying together as one word. I tried to separate the "postal" from the "name" by enlarging the word values, but it's only separate the label that have to words in it.

I have tried to  change the placement of the labels, but they all moved as one word.

How can I separate this two labels?

In this demo map example I need to separate the "ND" from the "North Dakota".


texts: Testing point in polygon overlay performance of Intersect, Spatial Join and Near in ArcGIS for Desktop?

texts:
In answers and comments to Select points within a polygon and update attributes I think it is being suggested that Spatial Join (Analysis) and/or Near (Analysis) may be faster than Intersect (Analysis) for performing a point in polygon overlay to transfer polygon attributes onto the points that intersect them.

This surprised me because I have always expected that Intersect might be a bit quicker than Spatial Join, and I think it would be quite a lot quicker than Near.  Unfortunately, I do not have an Advanced level license of ArcGIS for Desktop so I cannot performance test Near.

Does anyone have a reproducible test that indicates Near and/or Spatial Join being faster than Intersect?

I am going to provide as an answer my test on Intersect vs Spatial Join where I create a 10x10 fishnet of polygons and a 100x100 fishnet of points covering the same extent (0,0,10,10) and then perform a vanilla Intersect and a vanilla Spatial Join (INTERSECT) of the two.