texts: Symbology for traffic direction into and out of city

texts:
I'm using QGIS 3.4.7. I have one vector layer with road lines and a second with city polygons (and a third with city centroid points if that is better). I have already used symbologies with arrows to confirm the road directions match the direction of traffic in real life. My plan is to color the roads based on whether their traffic flows either into or out of the city, so I was wanting to use something like this in the field calculator after selecting my target city from either the polygon or point version of the city layer:

CASE
WHEN distance(end_point($geometry), targetcity) - distance(start_point($geometry), targetcity) > 0 THEN 'out'
WHEN distance(end_point($geometry), targetcity) - distance(start_point($geometry), targetcity) < 0 THEN 'in'
WHEN distance(end_point($geometry), targetcity) - distance(start_point($geometry), targetcity) == 0 THEN 'no_flow'
END

However, I haven't figured out how to use the field calculator while referencing another layer. I would ideally like a solution without a lot of Python coding because I'm fairly new to QGIS but I will learn it if I have to.


texts: Reclassify Slope using SAGA Raster Calculator

texts:
I am looking for the right formula for reclassifying a slope raster.
I need to switch the values 1 = bad, 10 = good to 1 = good, 10 = bad.
This Formula worked fine using the normal raster calculator:
10 - %"Recoded@1"% + 1
but I can't figure out the right formula using the SAGA or Gdal Raster Calculator for my processing model.
Any ideas!?


texts:
I am using the landsat BSL (Bare Soil Line) package from R https://rdrr.io/cran/landsat/man/BSL.html, and I am wondering why the code sample is using llimit 0.005 and not 0.01:

  BSL(band3, band4, method = "quantile", ulimit = 0.99, llimit = 0.005)

Source code is at https://rdrr.io/cran/landsat/src/R/BSL.R

I looked at the code and documentation, and documentation says: "The quantile method takes the lowest set of points, those with a NIR/red ratio less than the llimit-th quantile." Is there a reason why llimit parameter is set to 0.005 and not 0.01?  Does 0.005 set up a 2 tail region, one at 0.005 and one at 0.99?

I emailed the author of this package but I didn't get a response about this question so I'm hoping someone on GIS SE might know.


texts: Limit date range search when using Group Filter widget in Web AppBuilder Dev Edition 2.19

texts:
I am using the Group Filter widget with date fields. There is a point layer and a polygon layer in the web app. I don't want any restrictions on the point layer date range. However, the polygon layer has millions of polygons, so I want to limit the date range search to 30 days for the polygon layer. I also don't want the user to be able to make the polygon layer visible without first doing a filter (via the group filter widget).
I have set up the Group Filter widget like so (see screenshot). I am using WAB Dev Edition 2.19.
Currently, I have it set up so that when the user conducts a group filter, it automatically makes visible the point layer but not the polygon layer (code is in the Widget.js file). The user still has the option to turn on the polygon layer, but this is still causing lag issues on our system. So that is why I want to implement the above checks.
this._applyFilter(layer.layerObject, newExp, false);
//}
//layer.layerObject.setVisibility(true);
if(layer.id === "Poly_layer_Portal_5214"){
  console.log(layer)
  layer.layerObject.setVisibility(false);
}
else{
  layer.layerObject.setVisibility(true);
}


texts: Piping output of one processing algorithm as input into another algorithm in QGIS Processing script?

texts:
I created a processing script in QGIS 2.18 that connects two processing algorithms (first clip and than dissolve). The output vector from the clip is supposed to be an input into the dissolve.

The code of my processing script is provided here:

##input_layer=vector
##clip_layer=vector
##clipped_layer=output vector
##processed_layer=output vector

clip_area = processing.runalg("qgis:clip", input_layer, clip_layer,
clipped_layer)

fin_area = processing.runalg("qgis:dissolve", clip_area, "true", None,
processed_layer)

The screeshot of processing script UI when I try to run it is below:

Both of the algorithms work correctly standalone, but when I try to pipe them it does not work as intended. I only get clipped_layer as output and the following error notification in Python Console:

There are several questions adressing related issues like Pipe custom processing algorithm output into another in QGIS 3.0 but none really answers my question.

Does anybody know how to correctly pipe algorithms in the example processing script provided above?

I know it could be solved by using Processing Modeler but since I plan to expand the script further, please, don´t suggest this as solution.


texts: How to add a band with conditional value to a raster with PostGIS

texts:
I'm trying to convert a vector polygon into a raster, with pixel values that indicate the area of the intersection.  I have a query that gets the area of intersection between an empty raster and the polygon:

SELECT ST_Area(ST_Transform(geometry, 26986)) / 1000000 as area, x, y
FROM (
    SELECT ST_Intersection(ST_Transform(geom, 4326), wkb_geometry) as geometry,
    x, y
    FROM (
        SELECT (ST_PixelAsPolygons(rast)).*
        FROM empty_raster
    ) raster_pixels,
    vector_polygon
    WHERE region_name = 'Calgary-Buffalo'
    AND ST_Intersects(ST_Transform(geom,4326), wkb_geometry)
) intersects

Which gives me the output:

  | area  | x  |  y  |
1 | 4.998 | 79 | 203 |
2 | 1.535 | 80 | 203 |

The empty_raster table is an empty raster with width of 140 and height of 248.

What I need to do is to add a band that sets values to 0 except for the above x-y values.

I've thought to try something like this:

SELECT ST_AddBand(rast, '32BF'::text, area)
FROM empty_raster, see_above_query

But I don't know what my WHERE clause would be.

I have found that I should create the band with 0 values first:

UPDATE empty_raster
SET rast = ST_AddBand(rast, 1, '32BF'::text, 0)

And update the values with ST_SetValues or ST_SetValue.

Now, when I try to update the values, it only updates the first value:

UPDATE empty_raster
SET rast = ST_SetValue(rast, 1, x, y, area)
FROM see_first_query

Am I using this function incorrectly?  Do I need to do a for loop to generate the proper SQL?


texts: Resizing pixel resolution after interpolation using ArcGIS Desktop?

texts:
I did kriging spatial interpolation using climatic point file.

After doing interpolation, how can I resize the pixel resolution to pixel resolution of the NDVI data.

NDVI image's pixel resolution is 250m.


texts: QGIS Print composer Label Icon in Legend

texts:
In a legend of a map I want to label the icon with a number, just as it appears in the map. As example, I have different parcels and every parcel has a number. In the map, the parcels are labelled. In the legend, there only appears the border of the parcel, but not the label inside. Is it possible to get just a random number inside of the icon within the legend?

Hope you understand my question.


texts: How to publish SQL Server Spatial views to GeoServer?

texts: I have a web mapping service created on ESRI infrastructure (ArcGIS Server (WFS and WMS) with editable data from SQL Server Spatial. I am trying to switch from ArcGIS to Geoserver. I have published the layers from QGIS using GeoServer Explorer Plugin. But all the SQL Server Spatial tables are getting converted to shapefiles. The plugin allows to connect to PostGIS as in the linked page, but currently all my data is on SQL Server db and heavily related to quite a many attribute tables. Also since the web application has the editing workflow, so publishing static shapefile won't work for me. Is the way to publish views from SQL Server spatial to GeoServer is to publish views from PostGRE linked to SQL server? Are there any fallback for the approach - performance etc.? I am open to experiment with any other open source product if there is any setback for this approach.


texts: Add images and other text in QGIS 3.0.2 (Windows)

texts:
Today is my second day of using a GIS software. I have imported latitude, longitudes from csv into QGIS. I want that whenever users click on any location, they can see a pop with info: Country and State only, not all the fields that I see now.

Also, I want to add a field that will contain a hyperlink. In addition, I want to that users see the locations' image post clicking (or hover over). For this I am following this tutorial, but not helpful. Working file is here.


texts: How to force QGIS 3.6.1 create raster's pyramids (overviews) when I drag&drop them in QGIS?

texts:
There is a build-overviews-pyramids tool to create pyramids but I want to create them automatically when I add (drag & drop) raster layer into QGIS interface.

I have found parameters that is looking very similar to what I need but perhaps it only create pyramids when I use GDAL tools.


texts: Convert lat/long four corners data to raster

texts:
I have a data set in the following format.each row represents the coordinates of the four corners of a cell and the last column is the value associated to the cell.

lon1 lat1 lon2 lat2 lon3 lat3 lon4 lat4 TotHg
========================================================
-44.48, 49.88,  -44.69, 50.51,  -43.71, 50.64,  -43.5,  50.02, 8.50

-42.53  50.14   -42.71  50.77   -41.71  50.88   -41.54  50.25   9.02

-40.56  50.35   -40.71  50.98   -39.7   51.08   -39.56  50.44   9.21

-38.56  50.52   -38.68  51.16   -37.66  51.23   -37.56  50.59   8.75

-36.55  50.64   -36.64  51.28   -35.61  51.33   -35.54  50.69   8.53

I would like to convert it to a raster to be used in GIS. At the moment I do not know the coordinate system but I am investigating. All I know is that the resolution is supposed to be 2.5°. How can this be done in either R or in ArcMap? I have read a similar thread but I am not sure if it is the same situation since I have the coordinates for the four corners of each cell.this is the other question:
How can I convert data in the form of lat, lon, value into a raster file using R?

I have to resuscitate this question because I have received some information about projection. I am not sure if this will change the answer or help in any way. I am still quite confused at the moment. Here is the information I have received:

My model is define on a 2-d plane by using a polarstereographic projection with a griddistance of 150 km at 60N.
The projection is define by:

X_in_grid_units=x_coordinate_of_northpole+6370km/150km*(1+sin(60N))*tan(45-latitude/2)sin(longitude-(-32))
y_in_grid_units=y_coordinate_of_northpole-6370km/150km(1+sin(60N))*tan(45-latitude/2)*cos(longitude-(-32))

where 6370km is the used diameter for the projection, 150km as told before is the griddistance at 60N, -32=32W is central longitude for the projection, eg. the longitude parallel with y-axis.
The griddistance as a function of latitude is defined as d=150km*(1+sin(latitude)/(1+sin(60N)), eg. 160.77km at northpole, and 80.38km at equator.


texts:
Currently, I am using Google Maps API with HTML5 geolocation to create a POC about asset tracking via a web service. The user opens the web browser and we track it from there and display it on an admin map. I have been testing it out and the results are not great. A sample picture of my commute to work is below:

As can be seen, the dots are scattered all over the train tracks and not very accurate. I am trying to make it so that it is very accurate, within 10m of the actual position. Having done some research, HTML5 Geolocation is not that great, and I have read about using Apache Cordova/ Phonegap and native geolocation off the phone with google maps for more accuracy. Can someone verify that this is the case?

Eventually, we are planning on a native app, but for now, it is web based,but if I can make a strong enough pitch, we might go to a hybrid until native app dev is done.

In terms of what my setup is, I use navigator.getCurrentPostion() to get the initial position and then run watchPosition with an update to the DB every 30s. Both are set to high accuracy mode, but the results do not indicate this.

EDIT: took a look at Google's geolocation off of their site, and here are the results :


texts:
I have 3D buildings in .kmz file. I can see them in Google Earth, but I cannot import them to ArcGIS. What I've tried until now:

I saved it as a .kml file in GoogleEarth and then try to convert the .kml with DNR Garmin to a .shp file. This doesn't work for me. I get zero objects.
I also tried Kml2Shp Online but that also doesn't work. I uploaded my .kmz file successfully but then I cannot continue - it says there is no objects...

Any ideas how to solve this problem?


texts: Making heat map based on values assigned to polygon using ArcGIS Pro?

texts: I am very new to ArcGIS Pro and I have loaded a SHP file with polygons and associated 'scores'. I would like to create a heat map based on these values e.g. a value of 1 would colour the polygon red, a value of 5 would colour it orange and a value of 10 would colour it green.


texts: What are my options in terms of implementing a WPS server?

texts: I know about PyWPS, are there any good alternatives?


texts: Using kriging to extrapolate values outside of sampling polygon in ArcMap?

texts:
I have data points with values along 4 perpendicular transects. The kriging tool that I am currently using interpolates within the polygon created by the 4 transects, essentially a rectangle. I am looking for a way to extend the output of the kriging outside of the rectangle created by the 4 transects.

Is a way to do this by creating a linear interpolation for the area outside of the polygon created by the transects?


texts:
I am trying to set up Python so that I can execute ArcPy scripts outside of ArcMap. I'm new to Python, however I have previously been able to execute Pythons scripts from the cmd prompt (PowerShell), however I can't execute any ArcPy related scripts (maybe I'm misunderstanding).

First thing is that when I type "python" in the cmd, I am getting a different version to the one ArcMap 10.3 is using, which is 2.7.

I have set the path in Windows Environment Variable -

I'm using a textbook and seemed to have followed the instructions, so i'm not sure where I've gone wrong.


texts:
I´m working on a project where I have to draw historical shipping routes of trade, whaling, etc. ships from the XIX century.

I´ve done some research and checked Maritime route prediction avoiding land masses but I want to be sure that I am going in the right direction because I have to learn this tool from scratch:

To draw the routes: I have mainly the origin and the destination point and sometimes intermediate harbours. I need a way to approximately draw from data as currents, winds, coastline and landmasses a route of how the vessels got from one place to other. I am currently learning postgis and postgresql. For my goal, is Pgrouting the most useful and simple way? This question> Create a route by using points (pgRouting)
still unanswered makes me think is not.


texts: ArcGIS C# RubberPoint.TrackExisting not working

texts:
I have a simple marker symbol I have drawn on the map. I am trying to use a rubber band to drag-move it to a new location. I think I am close because when I use .TrackNew instead of .TrackExisting the object with move a little on the screen.

The simple marker symbol:

    public void DrawControlPoint(Point3D screenPoint, ISpatialReference mapSpatialReference)
    {
        IScreenDisplay screenDisplay = ((IActiveView)MapController.Map).ScreenDisplay;
        IGraphicsContainer graphicsContainer = GetGraphicsContainer();
        if (graphicsContainer != null)
        {
            IDisplayTransformation displayTransformation = screenDisplay.DisplayTransformation;
            IPoint mapPoint = displayTransformation.ToMapPoint(screenPoint.X, screenPoint.Y);
            mapPoint.Z = screenPoint.Z;
            mapPoint.SpatialReference = mapSpatialReference;

            ISimpleMarkerSymbol simpleMarkerSymbol = new SimpleMarkerSymbolClass();
            simpleMarkerSymbol.Color = EsriConversionUtility.ConvertToRgbColor(Colors.Blue);
            simpleMarkerSymbol.Outline = true;
            simpleMarkerSymbol.OutlineColor = EsriConversionUtility.ConvertToRgbColor(Colors.Red);
            simpleMarkerSymbol.Size = 9;
            simpleMarkerSymbol.Style = esriSimpleMarkerStyle.esriSMSDiamond;

            IMarkerElement markerElement = new MarkerElementClass();
            markerElement.Symbol = simpleMarkerSymbol;
            var element = (IElement)markerElement;
            element.Geometry = mapPoint;
            IElement elementToAdd = element;
            graphicsContainer.AddElement(elementToAdd, 0);
            MapController.PartialRefreshMap();
            if (ActiveMap.Instance.IsInitialized)
            {
                ActiveMap.Instance.Map.IsDirty = true;
            }
        }
    }

The drag move code:

    public void DragMoveControlPoint()
    {
        IElement element = _selectedElement;
        IScreenDisplay screenDisplay = ((IActiveView)MapController.Map).ScreenDisplay;
        IGraphicsContainer graphicsContainer = GetGraphicsContainer();
        if (graphicsContainer != null)
        {
            IPoint point = element.Geometry as IPoint;
            IGeometry geometry = point;

            screenDisplay.StartDrawing(screenDisplay.hDC, (Int16)esriScreenCache.esriNoScreenCache);
            IGeometry newGeometry = geometry;
            var rubberBand = new RubberPoint();
            //newGeometrym = rubberBand.TrackNew(screenDisplay, null);
            bool ok = rubberBand.TrackExisting(screenDisplay, null, newGeometry);
            screenDisplay.FinishDrawing();

            graphicsContainer.DeleteElement(element);
            element.Geometry = newGeometry;
            graphicsContainer.AddElement(element, 0);

            MapController.PartialRefreshMap();
            if (ActiveMap.Instance.IsInitialized)
            {
                ActiveMap.Instance.Map.IsDirty = true;
            }
        }
    }


texts:
I have inherited a QGIS system from a former employee and my predecessor built a custom plugin for it. I now need to remove a single checkbox option completely from the UI for this plugin. Is there a simple way to do this?

Thanks!


texts: Function that samples region for daily images within multiple years

texts:
I am trying to sample daily meteorological data from 1990 to 2018 within a region and export as a CSV file. How can I write a function within a function to filter first by year and then by day of year? Here is the syntax that I have so far. Is there is a faster way to do this?
var dataset=Gridmet.filter(ee.Filter.date('1990-01-01','2018-12-31'));
var bands = ['tmmx','rmax','rmin','tmmn','srad','th','pr'];

 var imgCol = ee.ImageCollection.fromImages(
                  ee.List.sequence(1991,2019).map(function(year){
 var data1 = dataset.filter(ee.Filter.calendarRange(year,year,'year'))
  return data1.filterBounds(CA).set('year',year);

  var imgColPerday = ee.ImageCollection.fromImages(
                  ee.List.sequence(1,366).map(function(day){
  var image = data1.filter(ee.Filter.calendarRange(day,day,'day_of_year'))
        .select(bands).mean();
  image=image.addBands(ee.Image.pixelLonLat());
  return image.set('day_of_year', day)}));// set the day
 }));

 var data =  ee.FeatureCollection(imgCol.map(function(image){

 // calculate data
  var MeteorCA= image.sampleRegions({
  collection: CA,
  scale: 4000
  }).map(function(feat){
  return feat.set('day_of_year',image.get('day_of_year').set('year',image.get('year')));
  });

   return MeteorCA;
   })).flatten();

   //Exporting data to google drive
   var polyOut = MeteorCA.select(['.*'],null,false);

   Export.table.toDrive({
    collection: polyOut,
    description:'MeteorCA',
    folder: 'MeteorCA90',
    });


texts: How to calculate elevation differences between points and terrain model?

texts:
I have a points layer with a Z (elevation) field and a terrain model (points layer or polygons layer -triangles- or eventually raster).
I'm looking for a way to select the entities of my points layer that are at least 10m under the terrain model.
Any way to do this in QGIS 1.8 ?


texts:
I want to visualize polygons in Leaflet. I do this via L.geoJson, no problem here.
But each polygon shall has it's own set of Markers, which should be dispayed, when the polygon is clicked and hidden, when a place outside of the polygon is shown.
Example: Marker 1, with lat1,lon1,description1 and Marker 2 lat2,lon2,description2 belong to Polygon A. Marker 3 and Marker 4 belong to Polygon B.
The default view is the map with the two polygons. No Markers. When I click on Polygon A, Marker 1 and Marker 2 should be displayed. When I click on Marker 1, description1 should be shown in a bubble or something, that comes out of Marker 1.
When I click on Polygon B, Marker 1 and Marker 2 should disappear and Marker 3 and Marker 4 should appear.
When I click somewhere on the map outside of the polygons, all Markers should disappear.
I created a featureGroup for the Markers and I can call clearLayers() to hide all Markers. But how do I add/create the markers so that each Marker is assigned to a polygon?
I can also add the Markers' information (lat, lon, text to display when clicked) as Features of the Polygon in geojson but then?
This is my code so far:
var map = L.map('map').setView([51.505, -0.09], 13);
var markerGroup = L.featureGroup().addTo(map);
var polygonGroup = L.featureGroup().addTo(map);

// clear all Markers when a place outside of a Polygon is clicked.
map.on('click', function(e) {
  markerGroup.clearLayers();
});

// Create the Markers of the clicked polygon
polygonGroup.on('click', function(e) {
  // this is needed, so that the onClick listener of map does not fire.
  L.DomEvent.stopPropagation(e);
  // Add Markers based on "properties" of a "feature" inside a "FeatureCollection" in a GeoJson File.
});

L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
  maxZoom: 19,
  attribution: '&copy; OpenStreetMap'
}).addTo(map);

var geojson = L.geoJson(segments).addTo(polygonGroup);

And this is the content of my geojson file:
var segments = {"type":"FeatureCollection","features":[{"type":"Feature","id":"30818214, 30437276.0","properties":{
"type":"Street",
"markers":[[[13.271,52.451],"text to display on this marker"],[[13.46,52.45],"text to display on this other marker"]]
},
"geometry":{"type":"Polygon","coordinates":[[[13.31650670833504,52.6062861338089],[13.316301432085778,52.60624094293158],[13.316165411632992,52.606064409617225],[13.316274301853358,52.60588107437412],[13.316564314460374,52.60579833204233],[13.316865565538109,52.60586465139],[13.31688028977927,52.60588376076443],[13.317085565006439,52.60592895087393],[13.317221588723905,52.60610548312916],[13.31711270191529,52.60628881921736],[13.316822687123965,52.60637156288567],[13.31652143261743,52.6063052434477],[13.31650670833504,52.6062861338089]
]]}},
{"type":"Feature","id":"30818215, 30437277.0","properties":{
"type":"House",
"markers":[[[13.271,52.45159],"text to display on this marker"],[[13.46,52.45],"text to display on this other marker"]]
},
"geometry":{"type":"Polygon","coordinates":[[[13.31650670833504,52.6062861338089],[13.316301432085778,52.60624094293158],[13.316165411632992,52.606064409617225],[13.316274301853358,52.60588107437412],[13.316564314460374,52.60579833204233],[13.316865565538109,52.60586465139],[13.31688028977927,52.60588376076443],[13.317085565006439,52.60592895087393],[13.317221588723905,52.60610548312916],[13.31711270191529,52.60628881921736],[13.316822687123965,52.60637156288567],[13.31652143261743,52.6063052434477],[13.31650670833504,52.6062861338089]
]]}}]};


texts: pgr_alphashape returns a self intersecting polygon

texts:
After using the function pgr_drivingdistance() I am attempting to make a polygon with pgr_alphashape. Following the example found in the documentation I should do the following: https://docs.pgrouting.org/2.4/en/pgr_alphaShape.html#pgr-alphashape

--Create a table with 1 Polygon where the agg_cost is lower then 5
DROP TABLE IF EXISTS temp.polygon;
CREATE TABLE schema.polygon as
SELECT ST_MakePolygon(ST_AddPoint(foo.openline,
    ST_StartPoint(foo.openline)))
    FROM (
SELECT ST_MakeLine(points ORDER BY id) AS openline
FROM (
    SELECT ST_MakePoint(x, y) AS points, row_number() over() AS id
    FROM pgr_alphaShape('
        SELECT
            source::integer as id,
            ST_X(the_geom)::float as x,
            ST_Y(the_geom)::float as y,
            agg_cost as cost
        from
            temp.exampleresult as kar
        INNER JOIN schema.table_vertices as vert
        ON kar.source = vert.id
        where agg_cost < 5;')
    ) AS a
) AS foo;

It works and it generates a polygon. However the polygon is self intersecting and has a weird shape. Where should I be looking to debug?

So the goal is to create after a nicely filled network with pgr_drivingdistance() a polygon which is based on the agg_cost.

This is what I am aiming for: https://anitagraser.com/2011/02/12/drive-time-isochrones/

I suppose the sequence of my points is not correct.

Ok so without my join of the agg_cost it works perfectly for all vertices. But with the join it does not:

DROP TABLE IF EXISTS temp.polygon;
CREATE TABLE schema.polygon as
SELECT ST_MakePolygon(ST_AddPoint(foo.openline,
    ST_StartPoint(foo.openline)))
    FROM (
SELECT ST_MakeLine(points ORDER BY id) AS openline
FROM (
    SELECT ST_MakePoint(x, y) AS points, row_number() over() AS id
    FROM pgr_alphaShape('
        SELECT
            source::integer as id,
            ST_X(the_geom)::float as x,
            ST_Y(the_geom)::float as y,
            agg_cost as cost
        from
            schema.table_vertices')
    ) AS a
) AS foo;


texts: How to use Mapserver+Openstreetmap for tracking service effectively?

texts:
I have got Asset tracking installations for many clients.
The availability of GIS map data for every client location seems to be the tough part. I have got a MapServer WMS installation.

Can I use OSM data behind mapserver? Also I need to timely update with latest OSM data.
Is there any mapfile (MapServer map style configuration file) for OSM data?
And is PostgreSQL/PostGIS also a better option to include?

Update: I followed exactly the instructions on http://trac.osgeo.org/mapserver/wiki/RenderingOsmDataUbuntu as advised by mapper. Everything worked fine. Got the osm-gooogle.map.
But, when I browsed http://localhost/cgi-bin/mapserv?map=/var/www/osm-demo/mapserver-utils-svn/osm-google.map&mode=browse&template=openlayers&layers=all
There was few things to change in the generated mapfile. Mainly, the IMAGEPATH & IMAGEURL were inside METADATA tag. Moved them outside to WEB tag and everything worked fine. It started drawing Google styled map from OSM data stored in local PostGres DB.


texts: Raster calculator and styling with QGIS Analytics?

texts:
I am working with a Radiometric TIFF and trying to find a way to make finding anomalies within the raster quicker.

The image is an orthorectified mosaic of many images of a concentrating solar site. The site is comprised of parabolic mirrors with a glass tube that runs down the middle. As tubes break or start to fail, they will become hotter than normally operating tubes.

When a tube is failed it is fairly obvious to the eye when looking through the orthomosaic. However, with 100,000 tubes over 2000 acres, I am looking for a way to have the raster styled or to clip out only the failed tubes based on their thermal fingerprint. A failed tube, based on the altitude flown, is approximately 10 pixels long, 1-2 pixels wide and are significantly hotter than the pixels surrounding it.

Issues however, are that these images are taken over multiple days and varying temperatures. 1/2 of the tubes are on the "return side" and are cooler and 1/2 are on the "heated" side and are hotter. Hot tubes found in the morning may be the same temperature as normal tubes in the afternoon. Surrounding ground structures are also sometimes in the same temperature range as the failed tubes. Even though temperatures vary throughout the site - the fingerprint is still the same. a connected string of 10 pixels or so hotter than the adjacent pixels.

Is there a way using the raster calculator to enter a formula that would match this fingerprint? Is there a way to use styling to color only those temperatures in the failed range and clip out all other?

I have a shapefile of the entire site that I could use for a mask with polygon features for every tube I can use as a constraint to the calculation if that would be possible?

I added two images... because I can't add more than that :(

1) "broken" is a broken tube showing the pixel value for the broken tube

2) "adj_normal" is a normal tube next to the broken tube

3) would have been an image of the lighter shade tube under the broken that is failed, but not broken. It's pixel value was 7714

4) would have been a general site image not zoomed in to show a larger section but really has no baring on the problem, was just for reference.


texts: How to export attribute table to Excel using ArcGIS Desktop?

texts: In previous versions I could just export as a .DBF file, open it with Excel and then convert it to an Excel file. However, now in ArcGIS 9.3.1 the .DBF file can no longer be open by Excel. I have also tried saving as a text file but to no avail.


texts: Error when trying to access variables in the GRASS startup program

texts:
I am able to create a new location using the Location Wizard, however I would like to create a GRASS location in a Python script, so I was figuring out this can be done from the console.

In the Grass startup webpage it says that this can be done with grass command. I tested grass -v, grass74 -v and grass76 -v, but message from the console is always the same "Der Befehl "grass" ist entweder falsch geschrieben oder konnte nicht gefunden werden" (the command "grass" is not correclty written or cannot be found).

I googled this and found that the problem might be due to spaces or special characters in the path, so I reinstalled GRASS in C:\GRASS_GIS_76. This did not solve the problem.

My OS system is Windows 10 and I installed GRASS GIS 7.6.1 (new stable).

What could I try next?


texts: What projection is this in GDAL classification

texts:
I need to re-project a PNG from a given projection to the Google maps projection.
I have a grid of the source projection, looks like a Mercator projection.
What is the exact name of this projection in the GDAL classification?

I need to re-project this exact region of the map, so if it is difficult to tell exactly, an any projection which would be a good approximation for this region would do.

Sorry for the dumb question, I'm a newbie.

Thanks!


texts: Geodatabase add current date and time to attribute

texts: Is it possible in Arcmap 10.0 to create a geodatabase feature class that has a date field which auto populates with the current date and time, every time a new feature is created? Is this possible without programming? Have checked ESRI help but nothing. The feature class is not in ArcSDE so it is not versioned.


texts:
I tried to follow these intructions on creating google-styled maps from osm-data. I have the SpatialLite data in QGIS (2.18) and now tried to load the styles provided by Anita. I open the properties window of each layer, click on "style" -> "load style" -> "load from file" and select the respective style. Irrespectively of the layer and style I choose, I get this error:

  unexpected character at line 35 column 91

I don't get this error when using other stylesheets. Nevertheless, I'd like to use Anitas styles, so I hope it is possible to figure out how to change her files in a way that make them usable with QGIS 2.18. However, I do not know how to edit these files.

Can anyone suggest a sensible starting point?


texts: Showing attributes of multiple parcels on one layout using DDP with dynamic text?

texts:
I am using ArcGIS 10.1 for making maps showing the distribution of parcels and their attributes in a given area. I am working with anywhere from 20 to 50 parcels in a set, sometimes more. I often use Data Driven Pages (DDP) and dynamic text to speed the process along. My problem is that sometimes there are multiple parcels in close proximity to each other. I would like to be able to show the attributes of multiple parcels on one layout using ddp/dynamic text. How can I do this?

Update: This is an example of what I most commonly work on. I used ddp and dynamic text for the attribute data on the right of the map. I wish to use dynamic text for attributes for all parcels in the layout. How can this be accomplished? Further explanations and examples can be provided if need be.


texts: Pixelated DSM using LAS Dataset to Raster?

texts:
I am using ArcGIS 10.2.2, and I am trying to create a DSM out of some LIDAR data. After checking the LAS Dataset properties, I see that I have an average point spacing of approx 0.3.

I tried creating a dataset with all the basic settings, just setting the sampling value to times three the point spacing, so I set 1.

The image seems ok, but when I add hillshading, I feel it seems kind of clunky. I am wondering if I use the dataset to raster tool right? Is there some approach that is "best practice"?


texts: Where do images that are exported thru "Project - Import/Export Export map to image get stored?


texts: Arcpy scripts freeze on reaching a certain level on big datasets

texts:
I have a geodatabase with 300 featureclasses with around 4.5 million features in each featureclass. I am updating each featureclass with a value calculated through an algorithm.

The IDLE ide for arcpy script freezes after running for 20 featureclasses. Then I need to forcibly close the script.

How can memory management be done for such a huge dataset?


texts:
My fresh QGIS 3.10.3 installation from  https://qgis.org/debian-ltr eoan main via apt install on Ubuntu 19.10 is missing (at least) one CRS. More specifically, I cannot find/select the EPSG: 3035 within QGIS. At the same time the console output shows an error message: ERROR 4: Unable to open EPSG support file gcs.csv.  Try setting the GDAL_DATA environment variable to point to the directory containing EPSG csv files.

Any suggestions what got messed up here or if that's a bug within QGIS/GDAL?

QGIS and OS versions

QGIS version | 3.10.3-A Coruña | QGIS code revision | 0e1f846438
-- | -- | -- | --
Compiled against Qt | 5.12.4 | Running against Qt | 5.12.4
Compiled against GDAL/OGR | 2.4.2 | Running against GDAL/OGR | 2.4.2
Compiled against GEOS | 3.7.2-CAPI-1.11.2 | Running against GEOS | 3.8.0-CAPI-1.13.1
Compiled against SQLite | 3.29.0 | Running against SQLite | 3.29.0
PostgreSQL Client Version | 11.7 (Ubuntu 11.7-0ubuntu0.19.10.1) | SpatiaLite Version | 4.3.0a
QWT Version | 6.1.4 | QScintilla2 Version | 2.10.4
Compiled against PROJ | 5.2.0 | Running against PROJ | Rel. 5.2.0, September 15th, 2018
OS Version | Ubuntu 19.10
Active python plugins | openlayers_plugin; processing; db_manager; MetaSearch

*UPDATE:

1) the environement variable GDAL_DATA is set to /usr/share/gdal; however this directory does not contain any gcs.csv file (it seems that gcs.csv files are not included any more in gdal: https://github.com/OSGeo/gdal/issues/1471)

2) The EPSG 3035 is also not found when using gdal directly:

gdalinfo --version
GDAL 2.4.2, released 2019/06/28

gdalsrsinfo EPSG:3035
ERROR 1: ERROR - failed to load SRS definition from EPSG:3035


texts: Join attributes by location (summary) returns NULL

texts:
I am trying to count points within polygons (311 calls by Census tract). Both shapefiles are EPSG:4326, and I've done define projection just to double-check.
Every time I run "Join attributes by location (summary)" all values come back NULL. If not the projection, what else can cause this, or, if it is the projection, what else can I do to check?


texts: Generating line offset, perpendicular using PostGIS?

texts:
I'm trying to generate two line offsets from a line but I'm having problems generating them as they should be.

Using PostGIS

select st_asgeojson(st_union(the_geom)) from (
select ST_Line_Substring(ST_OffsetCurve(line_geom,0.00005),0,1) as the_geom from table_name where id = 7382
  UNION
  select ST_Line_Substring(ST_OffsetCurve(line_geom,-0.00005),0,1) as the_geom from table_name where id = 7382
UNION
  select (line_geom) as the_geom from table_name where id = 7382) as a

the problem I'm seeing right now is the angle between the start/end of given line isn't perpendicular to the start/end of the other lines

Result:

Expected:

EDIT: lines' geometry and srid

"LINESTRING(-8.69371814905576 40.5676685897826,-8.69334901004008 40.567864446704)",4326
"LINESTRING(-8.6936947144568 40.5676244216994,-8.69332557544112 40.5678202786208)",4326
"LINESTRING(-8.69330214084216 40.5677761105376,-8.69367127985784 40.5675802536163)",4326


texts: Error updating firstChild.data when updating WebCapabilities

texts:
I'm trying to update the WebCapabilties to allow uploads for a service definition.  I've followed the instructions from ESRI's website, but when I try to set the keyValue.nextSibling.firstChild.data equal to "Uploads" I get the following error message...

AttributeError: 'NoneType' object has no attribute 'data'

doc = DOM.parse(sdDraft)
configProps = doc.getElementsByTagName('Info')[0]
propArray = configProps.firstChild
propSets = propArray.childNodes
for propSet in propSets:
    keyValues = propSet.childNodes
    for keyValue in keyValues:
        if keyValue.tagName == 'Key':
            if keyValue.firstChild.data == "WebCapabilities":
                keyValue.nextSibling.firstChild.data = "Uploads"

Inside the XML, I can see the what needs to be updated.

<Key>WebCapabilities</Key><Value xsi:type='xs:string'></Value>

Needs to be...

<Key>WebCapabilities</Key><Value xsi:type='xs:string'>Uploads</Value>

If I manually update the value, I can see the firstChild.data value, if I don't, then firstChild is set to None, and I can't access the .data property of the firstChild.


texts: Determining inflow outflow of polygon based on flow accumulation?

texts:
I'm working with ArcGIS 10.2 to try to determine if a given site has some sort of outflow for the water coming into the site. I'm using the ArcHydro tools and have calculated the flow direction and flow accumulation for every watershed in the state. In addition, I have a shapefile that has ~25 polygons (self-created) that depict each of my sites.

I have approximately four broad types of polygons: (1) Type I polygons have stream channels running through them and thus drain well. (2) Type II polygons have small-scale streams coming in (they don't show up because I've set my flow accumulation threshold > 250 pixels). However these polygons DO show a larger-scale stream flowing out of the polygon, and can also be considered to drain. (3) Type III polygons have no apparent streams coming into or out of them, and do not drain. Lastly, (4) Type IV polygons have streams draining in, but do no have streams draining out.

I'm trying to find a way to:

A) calculate the inflow vs outflow to show that Type I and II polygons have more water flowing out than they do flowing in, and thus drain. Versus Type III and IV polygons that have a greater inflow than outflow.
B) develop an inflow/outflow model that can then be used on any random polygon in the state to determine which type of polygon it would fall into.


texts: Exporting results of supervised classification using Ecognition

texts:
I am trying to export the result of the classification. I set the parameters as below. However, when I open the shapefile in QGIS, there's no category called "class name". Thus, I can't change the symbology.
How can I fix that, and present the classified map with QGIS?


texts:
I have Road Network which was selected based on some criteria. Below is image representing it. Each beginning and end of Road Segment is marked with cross:

Selection led to fragmentation of this network. How can I merge (ArcGIS, QGIS) all lines that form consistent continuity, like in image below (colors represent 3 lines which should result after merging). Result of merging – based on below example – should be one layer, with 3 rows, where each row represents grouped lines. Road Network has several attributes, but I am not interested in having them in final, merged output. My goal is to group connected lines:

I have tried few approaches e.g.: Using Dissolve and Multipart to Singlepart in QGIS? Join / merge lines that touch and intersect in QGIS and none of it worked.
Dissolve and Multipart-to-Singleparts (QGIS) is not working as expected when 3 lines (or more) are touching each other:


texts: QGIS mapserver - fails to display some raster images

texts:
QGIS 2.0.1 running on Linux, with mapserve support. Randomly, it seems, some raster images simply fail to display. No errors, just no image (well actually the mapserve returns a blank white image). I've tried with both ECW and TIF formats; some load, some don't (it's always the same ones that will or won't load).

The map files (and their associated images) load fine in the desktop version of QGIS.

The image files aren't particularly big - less then 30 MB for the ECWs, and under 100MB for the TIFs


texts:
I am trying to install QGIS 2.14 in OSX El Capitan using brew and osgeo/osgeo4mac repository. But when I tried to install it I got the following error

CMake Error at /usr/local/Cellar/cmake/3.6.1/share/cmake/Modules/FindQt4.cmake:1328 (message):
Found unsuitable Qt version "5.6.1" from /usr/local/opt/qt5/bin/qmake, this code requires Qt 4.x
Call Stack (most recent call first):
CMakeLists.txt:272 (FIND_PACKAGE)

Apparently there is a problem in version used for Qt 5, I have Qt4 also installed but qgis installation has dependencies of package qt5 so it installs it anyway.

What I tried so far:

Use --ignore-dependecy qt5 flag, but nothing happens, no installation and no error
Make a symbolic link through qmake files, but get the same error

I have recently installed QGIS in another OSX device with no problems


texts: Adjust center of the map based on marker position Openlayers

texts:
I have a map in which a moving marker is placed to represent a vehicle. I am using setCoordinate() method to move the marker.

I need to set the center of the map to the position of the vehicle marker whenever the marker reaches any one of the four boundaries of map. How do I know that the marker is at the boundary in order to set the center of the map and set the center whenever it is at the boundary.


texts: How to draw polygons from the python console?

texts: I've been visiting and revisiting the page on geometry handling in the PyQGIS Cookbook: http://documentation.qgis.org/2.0/en/docs/pyqgis_developer_cookbook/geometry.html but can't seem to figure out how to get the polygon to draw from the Python console. Can anyone help?


texts: Convert .dwg file into .dwf file using Teigha tool

texts:
I am using teigha tool to convert my dwg file into dwf file. Can anybody help to explain the procudure for the same.

Your prompt reply will be greatly appreciated!


texts: PostGIS Geometry Query returns "Error: Operation on mixed SRID geometries" only for certain values

texts:
I have PostGIS table with two geometry columns, both defined with SRID 4326. I can insert into the table without problem, using the following INSERT statement (where lng and lat are values passed in programmatically):

INSERT INTO pad_meta (
    uuid, created, updated, name, origin, radius, area, expiry, creator
) VALUES (
    $1, now(), now(), $2, ST_GeomFromText('POINT(lng, lat)', 4326), $3,
    ST_Buffer(ST_GeomFromText('POINT(lng, lat)', 4326), $4), $5, $6
)

But when i query for an intersection using ST_Intersects, depending on the value of the point I get ERROR: Operation on mixed SRID geometries.

For example, this query works:

SELECT * FROM pad_meta where ST_Intersects(
    'POINT(-122.334172173172 46.602634395263560)'::geometry, area::geometry
) ORDER BY created DESC;

And this errors out:

SELECT * FROM pad_meta where ST_Intersects(
    'POINT(-122.334172173172 47.602634395263560)'::geometry, area::geometry
) ORDER BY created DESC;

Note, they are identical queries except the value of the longitude. I have experimented with different values, but not identified a clear transition point between the queries that work and don't.

I think I am fundamentally misunderstanding something. For the moment, I have solved/corrected/worked around the problem by re-formatting the query to use ST_GeomFromText and explicitly specifying the SRID:

SELECT * FROM pad_meta where ST_Intersects(
    ST_GeomFromText('POINT(-122.334172173172 47.602634395263560)', 4326), area
) ORDER BY created DESC;

But I honestly don't really understand what the difference is, or if this truly "the" solution.

Why am I getting an error only for specific values, and what is the proper way to format this query?

Here is my table definition for reference:

CREATE TABLE IF NOT EXISTS pad_meta (
  uuid CHAR(32),
  created TIMESTAMP,
  updated TIMESTAMP,
  name VARCHAR(128),
  origin GEOMETRY(Point, 4326),
  radius INTEGER,
  area GEOMETRY(Polygon, 4326),
  expiry TIMESTAMP,
  creator CHAR(32),
  PRIMARY KEY (uuid)
);

I have also verified that there is only one type of SRID in the geometry_columns:

SELECT f_table_name, f_geometry_column, srid FROM geometry_columns;
f_table_name | f_geometry_column | srid
--------------+-------------------+------
 pad_meta     | origin            | 4326
 pad_meta     | area              | 4326

I have also seen Why does only one of the following queries throw an error?, but since I am already explicitly defining my geometry SRIDs when inserting into the table, it seems like that is not what is happening.


texts: Replacing all  values in the attribute tables of all feature classes in a FGDB with the blank (empty) ones

texts:
I need to iterate through a file GDB feature classes to replace  values found in all the columns of attribute tables with the blank ones (empty'').

Here is my code, that doesn't work. Truth be told I'm not familiar with programming languages, but couldn't implement this task in ModelBuilder, that is why it seems the only way.
What's wrong with it?

import arcpy

fc = "C:\Users\Vadim\Downloads\GML\MolokovGP_p10.gdb"
fieldList = arcpy.ListFields(fc, field_type="integer")

fieldObs = arcpy.ListFields (fc)
fieldNames = []
for field in fieldObs:
    fieldNames.append(field.name)
del fieldObs
fieldCount = len(fieldNames)

with arcpy.da.UpdateCursor(fc, fieldNames) as curU:
    for row in curU:
        rowU = row
        for field in range(fieldCount):
            if rowU[field] == None:
                rowU[field] = replacementValue

        curU.updateRow(rowU)

This error occurs:

Runtime error
Traceback (most recent call last):
  File "<string>", line 4, in <module>
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\__init__.py", line 1131, in ListFields
    return gp.listFields(dataset, wild_card, field_type)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\geoprocessing\_base.py", line 344, in listFields
    self._gp.ListFields(*gp_fixargs(args, True)))
RuntimeError: ERROR 999999: Error executing function.


texts: Removing duplicate labels within a given radius in QGIS

texts: In my map, every street segment is an entity. So each one is labelled as you can see in the picture. I'm trying to remove duplicate street labels within a given radius in QGIS, but I can't find how to do it. I know how to do it in ArcGIS, but in QGIS I cannot find the option.


texts:
I would like to create a map (of a part of the North Sea area of The Netherlands and Germany) with a few labeled markers at specific locations (lon/lat coordinates) in the North Sea. I have a current development version of QGis available; I think this will become 3.0. (I cannot use a released version; this is what I can install on my OS.) The purpose of the map is as an image in a scientific paper, to illustrate the location of some offshore structures relative to the surrounding countries. So the printed size will be 12 cm (hor) by 4 cm (vert), say.

I am a complete beginner using GIS software. I have tried to follow the tutorial at http://www.qgistutorials.com/en/docs/making_a_map.html, but the color was purple and I couldn't manage to set appropriate dimensions. I also had a look at How to create a printable map of labelled markers with OSM background?, but it didn't seem to apply to my QGis install and out-of-date in general. I furthermore tried to download some OSM data as detailed on https://wiki.openstreetmap.org/wiki/QGIS, but no data was actually downloaded because runtime error: Query ran out of memory in "bbox-query" at line 1. It would need at least 512 MB of RAM to continue..

So I am at a loss as how to continue. A task that I had assumed would be easy (‘just cut out a piece of an existing freely-licensed map, add some markers, done’) turns out to be extremely involved. It is not necessary that QGis is involved, but other GIS software is most likely not available to me. (I previously tried to use Marble for the job, but that didn't work either.)


texts: Using ModelBuilder iterator for calculating raster bands for multiple rasters in ArcMap

texts:
I have a raster for each year from 1985-2015, so in total 31, and I need to calculate the mean annual value for each of them. Every raster has 12 raster bands which represent the months. The calculation itself is easy, but I would like to iterate the process in the ModelBuilder. So far the model only has two parts, the iterator and the raster calculator. I've chosen the workspace for the iterator but i can not access the raster bands in the raster calculator, only the raster itself.
Is there a way I can choose/access the raster bands from the raster calculator in the ModelBuilder and iterate the whole thing?


texts:
I would like to give a shot to the Time manager in QGIS with some raster datasets.

I load 2 rasters (both single band, 32bit float).

I add the first one in the Time Manager window with the start date: 2008-01-01 / end date: 2008-01-31 , it s fine.
Then I add the second one with the start date: 2008-02-01 / end date: 2008-02-31.

After I enter the second one I have got the following error message:

Traceback (most recent call last):
File "C:/Users/b/.qgis2/python/plugins\timemanager\vectorlayerdialog.py", line 70, in add_layer_to_table
settings  = self.extract_settings()
File "C:/Users/b/.qgis2/python/plugins\timemanager\rasterlayerdialog.py", line 24, in extract_settings
return ls.getSettingsFromAddRasterLayersUI(self.dialog, self.tempLayerIndexToId)
File "C:/Users/b/.qgis2/python/plugins\timemanager\layer_settings.py", line 77, in getSettingsFromAddRasterLayersUI
result.layerId =  layerIndexToId[ui.comboBoxLayers.currentIndex()]
KeyError: 2

Anyone has already seen this error?

ps: I use QGIS 2.8.1


texts: NAME is not defined when defining class arcpy

texts:
I am having an issue I have not had before creating a class with a method as below;

class Points:

def __init__(self, inFC, buffDist, sector): #Must have a set of points and a value to buffer these by. Also, sector.

    self.inFC = inFC
    self.buffDist = buffDist
    self.sector = sector

def getCoords(self): #getting the coordinates of the points.

    fc = self.inFC
    fields = ['SHAPE@XY']

    coordsList = []

    with arcpy.da.SearchCursor(fc, fields) as rows:
        coordsList = [r[0] for r in rows]
    self.coordsList = coordsList
    del coordsList

The class is created and the method called from this code at the end of my script;

if __name__ == '__main__':

inFC = arcpy.GetParameterAsText(0)
buffDist = arcpy.GetParameterAsText(1)
sector = arcpy.GetParameterAsText(2)
outFC_data = arcpy.GetParameterAsText(3)

fcName = outFC_data.rpartition("\\")[2]
fcPath = outFC_data.rpartition("\\")[0]
outFC = arcpy.CreateFeatureclass_management(fcPath, fcName, "POLYGON")

pointobject = Points(inFC, buffDist, sector)
pointobject.getCoords()

This returns the error "fc is not defined" at the line fc = self.inFC.

OR

if I remove "fc = self.InFC" and replace fc with self.inFC within SearchCursor then I get an error at the line "with arcpy.da.SearchCursor(self.inFC, fields)" that states "self.inFC is not defined".

I have tried pasting the class definition directly into the python interpreter (ArcCatalog->Geoprocessing->Python) but the errors remain the same, so I do not even get to the point of creating an instance of the class, the error is within the actual code of my class and method definition.

What is the error in my code or my approach?


texts: Select by Attribute on Date field using wildcards?

texts:
I have a date field within a file geodatabase (10.1).

An example of a value in this field: 19/11/2012 08:49:56

I want to select all rows where the time is between 9pm to 11pm AND 6am to 9am.

I don't care what day, month, year it is.  So I want to make use of wildcards to ignore these values.

What would be the correct statement to put in the Select By Attributes to achieve this?


texts: Using input parameters to select layer by attribute

texts:
I am trying to create a Python-Toolbox with following Tool:

inputFeature: A User has to load a Feature
inputField: Then he has to to select a field from the Input Feature
inputValue: At the end he has to put in a value (no matter which datatype: date, string, numeric...)

The execute method now should select all rows which are equal to the inputValue.

I already had a working code, but just for numeric inputValues.
Another User recomended this code form a past post:

import arcpy

def buildWhereClause(table, field, value):
    """Constructs a SQL WHERE clause to select rows having the specified value
    within a given field and table."""

    # Add DBMS-specific field delimiters
    fieldDelimited = arcpy.AddFieldDelimiters(table, field)

    # Determine field type
    fieldType = arcpy.ListFields(table, field)[0].type

    # Add single-quotes for string field values
    if str(fieldType) == 'String':
        value = "'%s'" % value

    # Format WHERE clause
    whereClause = "%s = %s" % (fieldDelimited, value)
    return whereClause

if __name__ == "__main__":
    inputfc = r"C:\input.shp"
    outputfc = r"C:\output.shp"
    fieldname = "StudyID"
    fieldvalue = 101
    whereclause = buildWhereClause(inputfc, fieldname, fieldvalue)
    arcpy.Select_analysis(inputfc, outputfc, whereclause)

This is the link to the post: Including variable in where clause of arcpy.Select_analysis()
It sounds that this would be the perfect solution. But I have no idea how to implement the method in my Python-Toolbox and how to run it with Input-Parameters. Does somebody know how to put the code into the standard Python-Toolbox Codeblock? I mean shouldn't it be used in def execute?


texts: Mass reprojection of raster files in a loop [Linux]

texts:
I have a folder with numerous .asc raster files that each accompanied by a corresponding .prj file. The file basename is the same for both files. The projection however can vary between all of them (they are currently in various different UTM projection).
What I want to achieve is to standardize the projection of all those individual raster files.
My current approach is this bash-script

#!/bin/bash
echo Starting Conversion to WGS84
for f in *.asc;
 do
  echo Converting "$f";
  echo To "${f%.asc}_wgs84.tif";
  gdalwarp "$f" "${f%.asc}_wgs84.tif" -t_srs "+proj=longlat +ellps=WGS84";
  done
echo Conversion finished

However the script fails for each file returning:

  "Invalid request to pj_param, fatal"

So I would like to see a bunch of reprojected .tif (or any other format) files that all have the same projection.

While typing this I am also working on a python script, which however does not seem to be as straight-forward as bash.

call: gdalinfo grid2000.asc
It seems to fail even at the end?

Driver: AAIGrid/Arc/Info ASCII Grid
Files: grid2000.asc
       grid2000.prj
Size is 16, 16
Coordinate System is:
PROJCS["WGS_1984_UTM_Zone_30N",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-3],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (565146.201264380943030,680158.854243048001081)
Pixel Size = (30.000000000000000,-30.000000000000000)
invalid request to pj_param, fatal

call:
gdalsrsinfo xyz.prj # (Note that xyz.prj could be a different file)

PROJ.4 : '+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs '

OGC WKT :
PROJCS["WGS_1984_UTM_Zone_17S",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-81],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["Meter",1]]


texts: QGIS is flipping my rasters when I vectorize them

texts: I am using QGIS to vectorize some raster layers.  I was able to do this fine a week ago, but now every time I use the vectorize tool, the resulting vectors are flipped upside down.  Any suggestions?


texts:
I am working with the Modis version6 NDVI/EVI 16 day composite dataset (MOD13Q1). I am developing a automated algorithm to download, aggregate, and processes the dataset. Therefore, I do not want to use the MODIS specific tools, but use the gdal library to reproject the dataset to UTM 13N. However, I have found conflicting information regarding the correct projection information. spatialreference.org has several different projects:

http://spatialreference.org/ref/?search=MODIS

And here is the gdalsrsinfo for one of the HDF4 datasets:

gdalsrsinfo 'HDF4_EOS:EOS_GRID:"MOD13Q1.A2017049.h10v04.005.2017066041938.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days NDVI'

PROJ.4 : '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs '

OGC WKT :
PROJCS["unnamed",
GEOGCS["Unknown datum based upon the custom spheroid",
    DATUM["Not specified (based on custom spheroid)",
        SPHEROID["Custom spheroid",6371007.181,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]],
PROJECTION["Sinusoidal"],
PARAMETER["longitude_of_center",0],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["Meter",1]]

Can anyone provide some clarification to the correct SRID information?


texts: Vector layer invisible; attribute table empty but large shp and db files exist

texts:
Systems being used: QGIS 2.8.6. on Linux Ubuntu and also on Mac Sierra.

I have been sent a large shp file 13.4MB and dbf file 3.4MB plus the other 6 files for the shapefile vector layer (all probably created in ArcGIS). All these files are held together in the same folder and their names are Hcc_cropmarks.shp etc.

This data represents archaeological features seen in aerial photographs and the shape file should contain a very large number of 'line' features varying in size from tens of meters to several kilometers long. The owner of this data is not apparently having any problem seeing it themselves in their ArcGIS system and previous shp files they have supplied worked perfectly for me.

I used 'add vector layer' to add the shp layer to an existing project. It is the topmost layer and is invisible even if all the other layers are turned off. I used CRS EPSG:27700 OSGB1936/British National Grid) which is the same as my project.

I believe that this is the correct CRS/projection as opening the prj file with gedit shows the following:

PROJCS["OSGB 1936 / British National Grid",GEOGCS["OSGB 1936",DATUM["OSGB_1936",SPHEROID["Airy 1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],AUTHORITY["EPSG","6277"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4277"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","27700"]]

My problem is that the layer is invisible. 'Zoom to layer' causes a zoom out for the other layers in my working area (but they are still visible) but the cropmark layer is invisible.  I can open the dbf with Libre Office on the Mac (but not on the Linux box which gives error message "General error. General input/output error.). The data in the dbf is arranged as a table of 3 columns (headers are OBJECTID,N,10,0; CROPMAR_ID,N,10,0; Shape_STLe,N,19,11) and about 100,000 rows, sequentially numbered in the first column.

When I open the attribute table for the layer it shows the 3 headers (without the N,10,0 etc) but no other data - the table is otherwise apparently empty.

The layer does not appear to be making the connection between shp and dbf.

I need to know if the visibility problem is of my making or whether the data I have been sent is incomplete/corrupted, in which case I could go back to the creator and ask for a fresh copy - but obviously  I don't want to trouble them if it is my error.

Any suggestions for making the layer visible?


texts:
I am adding property names to a property boundary shapefile using QGIS 2.8.2.

I don't know what I have done to get several "qgis rocks" text boxes on my screen.

How do I remove them?


texts:
I don't really get the behaviour of QGIS new reports (3.X) ...

I inserted an attribute table in a report. this table is spanning over several pages. I would like to create a page footer with the page numbering (current page and total amount of pages).

That seems quite trivial ... but ... can QGIS reports do that ?

I'm not in the case of an atlas.


texts: Extruding building heights from 3D terrain model creates sloped buildings

texts:
How can i prevent sloped buildings when I extrude building footprints over sloped terrain? I am using FME to drape 2D features over terrain (SurfaceDraper), then using the Extruder transformer to create 3D buildings from OS Building Heights. Unfortunately the height value seems to apply itself to each vertex, which creates sloped building tops (roof geom is not relevant here) - example below:

Is there an extrusion method that uses a single vertex, or is there another method entirely?


texts:
After a great deal of research here is what is happening:
We select two points on a hydrography map with an island inbetween.  It is expected that the CostDistanceWithPolyline would return a line that traverses around the island.  We are now finding that the Cost Array that is computed ignores the existance of SOME of our NoData areas, but not all of them!  What is wrong here???  We cannot find the issue!!!

The hydrography map with island in question: note the value is NoData for a point I chose inside the island geometry.

The resulting shortest path.  Kinda hard for a boat to do that!

The cost array and it's value at the island location.

The call is so stupid for getting a cost array that I cannot determine how the call could be messed up.  As per code that ESRI did for us: we pass the hydrography map, the target point, and call GetCostDistanceFull

Here is our call to get the cost distance and backlink arrays.

    private static IGeoDataset[] ComputeFullRaster(IFeatureClass origin, RasterDataset costArray, IGPMessages messages)
    {
        List<IGeoDataset> results = new List<IGeoDataset>();
        IDistanceOp distanceOp = new RasterDistanceOpClass();

        // Create the cost distance
        messages.AddMessage("Creating the distance cost raster...");
        object objectMissing = Type.Missing;

        results.Add( distanceOp.CostDistanceFull(
            (IGeoDataset)origin, (IGeoDataset)costArray,
            true, false,
            false, ref objectMissing, ref objectMissing) );
        messages.AddMessage("Creating the backlink cost raster...");
        results.Add(distanceOp.CostDistanceFull(
            (IGeoDataset)origin, (IGeoDataset)costArray,
            false, true,
            false, ref objectMissing, ref objectMissing));
        return results.ToArray();
    }


texts: Distance between coordinates computation speed: conversion to UTM vs Great Circle

texts:
For a program I am writing I need to calculate the distances between coordinates. Seeing as I need to perform a large number of these calculations for each coordinate (thousands) it makes sense to convert to UTM, then calculate the distance (seeing as I am just comparing distances, I can even work with the square of the distance between the points).

At what point does this trade off make sense? As far as I can tell by reading the relevant calculations Greater Circle would be the faster algorithm, but obviously once you do the conversion to UTM you don't need to do it again. How many times do you need to run the Greater Circle algorithm before it makes sense to convert to UTM instead? (Assume that the accuracy isn't a big deal in my situation).


texts: Difference between DEM raster in QGIS 2.18 and 3.4.4

texts:
I came across a really strange problem.

Untill today i always worked in QGIS 2.18. Now i wanted to upgraded to 3.4.4 ( it was already downloaded but hadn't used it yet)

When opening my QGIS project everything looked fine except the styling of the DEM raster was wrong. I thought it was a small issue because the project was build in an older version and i just decided to reload the raster layer and restyle it.

After reloading the DEM i noticed that the elevation was wrong.

This is the raster in QGIS 2.18 :

And now in QGIS 3.4.4 :

And i really don't know what the problem can be.. I tried googling it but as far as i know i'm the only one with this problem.

The DEM is originally a .pts file converted to a DEM raster with autocad C3D.

Does anybody have a clue what might be the problem?

Edit:

QGIS 2.18 min / max value:

QGIS 3.4.4min/max value:

UPDATE:

Still no solution for my problem. Today i made a new DEM raster using Autocad Civil 3D. Min = MAX?

This are my raster settings

This is how it looks in the canvas :

This is the same DEM is autocad map :


texts: Getting QGIS Model Processing Raster Calculator Error?

texts: I am using "raster calculator" tool of the saga in my model. The tool correctly executes raster operation like "a * 2" in value formula option but the same expression is not executed pre-calculated value formula option.


texts:
I am using QGIS 2.4 on Windows 8.1. I have a project in which I have been digitizing from OpenMaps. Last night I was doing some more digitizing, saved, and closed the program. When I opened it this morning, my layers (.shp) were not projected properly, even though they were when I saved and exited. I closed, and re-opened, and they were projected properly (so weird thing #1). However, when I went to go back to digitizing, I noticed that I have some features of a layer showing up that do not exist anymore (I deleted them yesterday). These features show up only at certain zoom levels, and are not selectable. Screenshots attached.

Any idea what this means?

"Phantom" features, in the "Brush" layer. The first image is what it SHOULD look like, and the second image is what I see when I zoom out. This is happening in other layers as well.


texts: Calculating countries flown over on flight path

texts:
Say I have two airports given.

I can then get their coordinates.

Now that I have the beginning and end coordinates I can draw flight-path from the two locations.

Now my question is that given two end points is there a "preferred" way to calculate which nations the flight path flies over?

My two possibilities:

Approach 1: You "ping" the flight path at increments to see where the flight currently is using google's Geocoder. However I am under the impression this would cost a large amount of data.

Approach 2: Find some sort of polygonal data to represent the nations and see when countries' borders the flightpath intersect.

Now I'd prefer approach two; but I am not sure where to find good resources/data to use.

Any suggestions on a best approach?

Thanks guys.

EDIT 1:

So far information I've found are as follows;

Link to http://www.naturalearthdata.com/features/
Link below from the first responder
http://thematicmapping.org/downloads/world_borders.php

I notice that world_borders uses shapefiles; is there a way to read/import those into a java environment?

Any open-source engines that can read the data? I am assuming that the shapefile is what holds the polygonal country data?

Edit 2
For future readers:
I found that for programming Java; the link in the comment of the answer below had good information.

namely using GeoTools to complete the project they have a step-by-step guide for Netbeans as well as Eclipse that I will report back with if it solved my problem.

Furthermore, I am using their step-by-step guide in conjuction with data from NaturalEarth.


texts:
I want to generate an Atlas in QGIS using a fixed scale of 1:750 on each page. I am using a point file to define page location centers.
Both 'Margin around feature' and 'Predefined scale (best fit)' are greyed out but using 'Fixed scale' even when I set the scale to 1:750, when I switch pages the scale reverts to 1:3412.
I have added 1:750 scale in Project properties but nothing seems to change.
I am using QGIS 2.18.12


texts:
Using the initGRASS command in RStudio on Windows, I get the following error after I installed a new GRASS version on a new path:
Error : XML content does not seem to be XML: 'proj_get_authorities_from_database: Cannot find proj.db'
In addition: Warning message:
In system(cmd0, intern = TRUE) :
  running command 'g.proj.exe --interface-description' had status 5
Error in parseGRASS(cmd, legacyExec = legacyExec) : g.proj not parsed

Code snippet:
initGRASS("C:/GRASS GIS 7.8",
          home = tempdir(),
          mapset = "PERMANENT",
          override = TRUE)

For some reason RStudio computes something when I execute a GRASS function, but doesn't readRAST any result but showing the error message above.
The environmental variables, which I checked in the Grass shell, seem to be okay to me (considering: https://github.com/OSGeo/grass/blob/master/mswindows/env.bat).
Can anybody help solving the problem or at least explain the error message?
I just tried to use the Grass version 7.6 and it works. In the ongoing progress I will see, if using this version will last for my work. Until then I'm looking foward to your answers and comments.
Thanks
AJ


texts:
I see arcpy.ddd referenced in 3D Analyst code, but I don't see it documented anywhere.  What is this module and how does it relate to 3D Analyst?
For example:

Stack Profile

# Execute StackProfile
arcpy.ddd.StackProfile(profileLine, profileTargets, profileTable, graphName)
# Execute SaveGraph
arcpy.management.SaveGraph(graphName, outGraph)


texts: "Mosaic to new Raster" just stops working

texts:
I have been using the command "Mosaic to new raster" so far successfully. Now, without changing any option from before, it just doesn't work anymore. It is something it already did: pretending to work, but actually gets stuck until I close it and open it again.

It doesn't fail, it is just still showing it is working, but it has already done after 1 minute. This time it didn't not work like that any longer. I closed the software and it just give me back standard "error".


texts: Using a base layer from ESRI and converting to a regular layer and "cropping it" also using it for arcScene

texts:
How do I modify a baselayer from Esri?

I found a great one for my region (much better resolution than the one I was given) but it is way too big and I want to cut out the ocean part of the raster so I just have the city map.  It is sort of like what you would see with google earth. I would also like to use it for my ArcScene to make nice 3D images but it does not seem that is an option.

I am a bit of a beginner.


texts: How to create a map from the highest and lowest average values

texts:
Anybody can help me, how to generate the map from assumption data.
I have 100 points data (from GPS reading) showing the location of the bee nest. Unfortunately there is no record on how much volume they harvested from each nest. However, from my respondent, they gave me the two values from the highest and lowest average of the volume production, ex. 5.99 liter (highest) and 1.03 liter (lowest).
I want to make a map to see where is the highest possible production volume per unit area and the other showing the lowest possible production volume per unit area.


texts: "Geometry cannot have Z values" exception in ArcGIS Engine

texts:
I'm trying to put a list of point features (selected points) into a new FeatureClass. but when I add them to my FeatureClass I get the "Geometry cannot have Z values." exception on the line "buff.Shape = point;" .
can any one help?

        List<IFeature> features = GetSelectedFeatures(axMapControl1.get_Layer(0) as IFeatureLayer);
        IPointArray pa = new PointArrayClass();

        List<IPoint> p = new List<IPoint>();
        for (int i = 0; i < features.Count; i++)
        {
            p.Add(features[i].ShapeCopy as IPoint);
            pa.Add(p[i]);
        }

        IPoint point = new PointClass();
        point = pa.get_Element(0);

        point.SpatialReference = p[0].SpatialReference;

        IWorkspaceFactory workspaceFactory = new ShapefileWorkspaceFactoryClass();
        IWorkspace workspace = workspaceFactory.OpenFromFile(@"C:\", 0);
        IFeatureWorkspace featureWorkspace = workspace as IFeatureWorkspace;

        IFeatureClass FC = CreateFeatureClass("FC", featureWorkspace, pa.get_Element(0).SpatialReference,esriGeometryType.esriGeometryPoint);

        IFeatureCursor fCur = FC.Insert(false);
        IFeatureBuffer buff = FC.CreateFeatureBuffer();

        buff.Shape = point;
        fCur.InsertFeature(buff);

        fCur.Flush();
        System.Runtime.InteropServices.Marshal.FinalReleaseComObject(fCur);

        IFeatureLayer fl = new ESRI.ArcGIS.Carto.FeatureLayer();
        fl.FeatureClass = FC;
        var layer = fl as ILayer;

        axMapControl1.AddLayer(layer);

        axMapControl1.Refresh();


texts: Change marker color based on GeoJSON filter in Leaftlet

texts:
This is what I would like to do:

Display all features of a GeoJSON file using Leaflet. Provide a text box where a property can be entered, if any of the GeoJSON properties match with what is entered in the text box the marker icon color is changed.

Does anyone know of an example of this or an idea of how this should be done?


texts: Map algebra expression with dynamic raster layers

texts:
Attached is the model that i am working on.

The "To Value" and "From Value" are model parameters and hence subject to change each time according to requirements,
This leads to varying number of files generated each time. This variation in the number and names of output files in a workspace is creating the problem in hardcoding the following calculator expression in a subsequent model directly.

This is the expression in the second map algebra tool that i want to use. The input rasters are those generated from the model shown above.
Raster Calculator (output1+(output2/Output1)+(output3/Output2)+(output4/Output3)+(output5/Output4)+....)

I saw some effort to a similar question here,  Raster Calculator dynamic map algebra expression, but I didn't find it conclusive. (I didn't understand the advice to use the FOR loop there)

Help and advice please!
I am using ArcGIS-10.1


texts:
I have a project manager who really wants "chemboxes", basically labels/callouts with tabular data.  I've been using an older vbscript labeling expression, but they want something that is more aesthetically pleasing. I also require that it be either feature linked or easily updated/replicated with new data, but I don't mind some manual arrangement for clearer labeling/layout (I just don't have time to create individual tables).
The only thing I've really found that might match their requirements is geoCallout (http://www.superiorinfotech.com/products.html) but I can't seem to find any reviews.
Does anyone have any experience with the geoCallout?


texts: Fixing local orthorectification mismatches leading to small georeferencing errors?

texts:
My project uses good quality high resolution orthorectified imagery (from the government) as its reference base layer. I complement it for vegetation analysis with additional imagery from a different source, also orthorectified.
The georeferencing is pretty good overall, but there are local "registration differences" between the two in some places. They probably arise from different orthorectification algorithms and assumed DEMs used by the 2 vendors. See example below; the polygons are digitized pine trees from the reference leaf-off imagery that I consider "ground truth". The overlaid leaf-on imagery from the other source is misaligned by about 4 metres in this location.

What's the best approach to fix this, and workflow to ideally do it dynamically, in QGIS or compatible tools (e.g. GDAL, grass, other plugins)?
By "dynamically" I mean the following. I can of course re-georeference the supplemental imagery throughout my broad ROI, using a large number of GCPs and the thin plate spline method (which operates by local deformations). But this seems to be intended as a one-off, where I do it well once, and then use the modified (deformed) re-georeferenced image as my layer. It can be re-done at some point in the future with modified GCPs, but it's a hassle and time-consuming when the overall ROI/image is large.
What I'd love is a more "on the fly" version, where I can zoom into an area that becomes of interest (like the pic attached), fix up the errors there when I need to, and have QGIS (or whatever) "locally deform on the fly" just like it reprojects CRSes on the fly. I might later do the same in some other area of the image. I expect thin plate spline would likely be the algorithm under the hood, but I'm open to anything that works (piecewise affine likely good enough).


texts: GeoServer WCS 2.0.1 subset request data with different file dimensions

texts:
I'm trying get two sets of data using WCS 2.0.1 request from my GeoServer and do raster calculation using GDAL's gdal_calc.py. Both datasets in GeoServer are converted to the same coordinates and same resolutions, and the subset Long and Lat data are exactly the same in my two WCS requests.

How come I'm getting file with different dimensions?

The interesting thing is they are not off big, only the Y value is a unit bigger. The gdal_calc.py only takes raster layers with same dimensions, that's why I need to use same file dimensions.  Can anyone help me?

Below are my WCS requests, value of Long and Lat are posted by my front HTML page using JavaScript, since I was able to download data, the values are fine.

$url1 = 'http://127.0.0.1:8080/geoserver/ows?service=WCS&version=2.0.1&request=GetCoverage&CoverageId=Canada:canada_clay_250_sl1&subset=Long('.$minX.','.$maxX.')&subset=Lat('.$minY.','.$maxY.')';

$url2 = 'http://127.0.0.1:8080/geoserver/ows?service=WCS&version=2.0.1&request=GetCoverage&CoverageId=manitoba:crop_risk&subset=Long('.$minX.','.$maxX.')&subset=Lat('.$minY.','.$maxY.')';

Also, here are two screen captures of my published data info from GeoServer:

Data one for Canada_Clay

Data two for Manitoba Cropping History


texts: How to export a map into a low resolution image in QGIS?

texts:
The only way I know how to export a map view into an image is to resize the window and from there export it. But I want a separate ocean layer included in the image, so I click Zoom Full before making the ocean layer visible. I tried increasing the width of the panels to hide some areas I don't want included in the export, but the software stubbornly automatically adjusts the map zoom level, map extent, and the panel widths.

So I tried using Print Composer, which I try to avoid because of its complexity.

The only thing I can find is Add Map, but it must be dragged so that it is loaded into the canvas. Upon loading the map, the map extent and zoom level was very different from the QGIS main window.

Why can't I just paste it directly into the canvas so that the map extent is not modified?

I have several maps I plan to have the same resolution, but I could not find any plugins that would make exporting easier.

How can I export into a fixed map view into a low resolution image?


texts: How to cluster points based on a kernel density?

texts:
I have a large dataset with 36k points representing commercial land uses, each with a field containing the square footage. I have run a kernel density analysis on this dataset, producing a raster showing the density of commercial square footage over the entire metro area. I need to divide this raster into regions corresponding to local maxima, which I call a "center". I have already determined the locations of the centers, and now I need to do one of two things:

use a point clustering tool, such as "partitioning around medoids", to group the points into clusters around the centers I've identified. The problem with this method is that it is computationally intense, and even more so if I try to use a dissimilarity matrix to weight the points by size.
somehow divide the kernel density raster (which roughly resembles a terrain raster) into individual "hills" around each center. But I can't think of any tool for doing this.

This problem has plagued me for a while, and I hoped that I'd be able to perform the clustering method in R, but it is time consuming and I'm running out of time. Does anyone know of a simple method for either dividing density rasters into neighborhoods of intensity or for quickly clustering large datasets?


texts: How to make a WMS GetFeatureInfo request in Leaflet?

texts: I am trying to make a WMS GetFeatureInfo request using Leaflet. I can not find something in the documentation to describe this process. Can you please provide with an example demostraing how to do something like this?


texts: Creating and updating a field from calculation

texts:
I'm trying to create a field in a vector layer and then update the values based on the sum of two fields which are already present in the table, using the Python console (QGIS3.4). Also I'm still learning Python...
I originally started to develop this under 2.18, but I've recently had to move to QGIS 3 (which is part of the problem)- In 2.18 the first part of the script used to work  i.e. creation of a new field in the table.

The original code:

from qgis.core import *
from qgis.utils import *
from qgis.gui import *
from PyQt4.QtCore import QVariant
import processing
import os

layer_P = QgsProject.instance().mapLayersByName("BXO")[0]

SUM = QgsField('t_branch', QVariant.Int)
layer_P.addAttribute (SUM)
idx = layer_P.fieldNameIndex('t_branch')

layer_P.startEditing()

e = QgsExpression ('CB_Num_Ent' + 'CB_Num_Ext')
e.prepare( 't_branch'() )

for f in layer_P.getFeatures():
    f[idx] = e.evaluate( f )
    layer_P.updateFeature( f )

layer_P.updateField()
layer_P.commitChanges()

I've been looking at the API and seen that fieldNameIndex has been replaced with fields() .lookupField() or fields() .indexFromName() and since I've tried a couple of variations - for example:

layer_P = QgsProject.instance().mapLayersByName("BXO")[0]

SUM = QgsField('t_branch', QVariant.Int)
layer_P.addAttribute (SUM)

idx = layer_P.lookupField('t_branch')

layer_P.startEditing()

e = QgsExpression ('CB_Num_Ent' + 'CB_Num_Ext')
e.prepare( layer_P.fields() )

for f in layer_P.getFeatures():
    f[idx] = e.evaluate( f )
    layer_P.updateFeature( f )

layer_P.updateField()
layer_P.commitChanges()

but for the minute I'm still get the same type of error message
...AttributeError: 'QgsVectorLayer' object has no attribute 'lookupField'

So, firstly I'm trying to find out how to modify the script following the removal of indexFromName. And if that works, hopefully update my new field based on the calculation here :

> e = QgsExpression ('CB_Num_Ent' + 'CB_Num_Ext')

For info the script was adapted from the following :
Adding field and calculating expression with PyQGIS

PyQGIS: Multiply Fields and populate new field

How to fill fields with layer name in pyqgis


texts:
There was a table manager plugin in earlier versions I can't see anymore.

Is there a way to achieve this in QGIS 3.6?


texts: Creating features from JSON without writing to disk first using ArcPy?

texts:
In ArcPy, I want to use JSONToFeatures_conversion, but I'm working with data that I pulled from a web service. I have no need to save the data to disk. So, I thought I could do something like this:

fd = StringIO.StringIO(<my JSON string>)
arcpy.JSONToFeatures_conversion(fd, 'in_memory\\output')

arcpy definitely doesn't like that though. I'm considering using the tempfile module, but that's not an ideal solution for this. Temp files aren't really a good solution for this because of a few factors, from environmental constraints (permissions in some of the environments where it may run) to the size of the datasets (IO overhead).

doing this:

fs = arcpy.FeatureSet()
fs.load(<url to service>)

Isn't exactly ideal either because the dataset is composed of data from a couple sources, with different attributes that have to be reconciled.

Ultimately, the product of this operation will be passed in to arcpy.SelectLayerByLocation_management as the select_features parameter

Environment: ArcGIS 10.x (multiple machines running slightly different versions, but none of it's ArcPro, hence all Python2 not Python3)


texts: Deleting field in shapefile based on values in field using ArcPy

texts:
I'm trying to create a Python script to automate simplifying some feature classes (shapefiles).  The data table for my shapefiles look like this (all of the fields except the first four are SHORT):

FID    Shape*    GEOID10            EMS_ORI    A_CO01    A_CO05    A_CO06    A_CO11
___________________________________________________________________________________
0      Polygon   421110218001027    E56001     99        23        99        17
1      Polygon   421110218001231    E56001     99        24        30        15
2      Polygon   421110218001028    E56001     99        24        99        16
3      Polygon   421110218001035    E56001     99        24        99        99

and so on.  I am trying to find an easy way to delete an entire field if it only contains the value "99" for every row.  In this example, I would want to delete the A_CO01 field because all features in this shape file have a value of 99 there.  I do not want to delete A_CO06 or A_CO11, because at least one feature has a value other than 99 in each of those fields.  My shapefiles have thousands of features, so I tried to script this with Python in an instant Python window in ArcCatalog.  I'm only a novice at Python, so I know I'm missing something easy, but here's what I tried:
import arcpy

fc = 'E:/EMS500_Polygons.shp'
fields = arcpy.ListFields(fc)

for field in fields:
    if (field.name == "FID" or field.name == "Shape" or field.name == "GEOID10" or field.name == "EMS_ORI"):
        continue
    else:
        if (field.type = "Integer" and min(field) == 99):
            arcpy.DeleteField_management(fc, field)

This gives me a syntax error on line 10: if (field.type = "Integer" and min(field) == 99):.  What I'm trying to do is delete any field where the MIN value of the entire field is 99 (this value is a "placeholder"; valid values for this field range from 1 to 93).


texts: Is there a way to prepare a project on QGIS with many layers and to publish it under WMS service on GeoServer?


texts: How to restore a scarp surface grid to calculate material loss

texts:
I would like to calculate the volume of material lost from a scarp like the one shown in the following image:

Is there a way to do these calculations with open-source software such as GMT, QGIS and GRASS?


texts: I have forgotten my GeoServer login id & password. Is there any way to restore default password? Or is there any way to view the password?


texts:
Intent is to have two variables (QgsVectorLayer and feature.id() within said layer) as input for the Identify feature map tool in Qgis (2.6).
Wondering if a feature can be selected programmatically (via PyQGIS) to be highlighted in mapCanvas()?

I have reached 'writers' block after the following:

    layerID,featureID = [i for i in layerAndFeatureIDs.split("_")]
    layerRef = ftools_utils.getMapLayerByName(unicode(layerID))
    infoTool = QgsMapToolIdentifyFeature(self.iface.mapCanvas(),layerRef)
    ...?

Runs without error - but hopeful to find out if calling a specific feature within a designated layer is possible.
Either by way of a feature.id(), or even the centroid of said feature as a QgsPoint reference.

This command might be relevant.  It activates the tool in question, but is it possible after triggering to subsequently highlight or info-select programmatically?

        self.iface.actionIdentify().trigger()

I wondered if an alternative method would be to source the feature in questions centroid.
Yet upon deriving such - and (futily?) implacing that location into a simulated click event on mapCanvas(); no results were had.  Although, a print of 'working' did return...

    for feature in layerRef.getFeatures():
        if feature.id() != int(featureID):
           pass
        elif feature.id() == int(featureID):
           x,y = feature.geometry().centroid().asPoint()[0],feature.geometry().centroid().asPoint()[1]
    self.iface.setActiveLayer(layerRef)
    self.iface.actionIdentify().trigger()
    virtualMouseClick = QtGui.QMouseEvent(QtCore.QEvent.MouseButtonPress,QtCore.QPoint(x,y), \
                        QtCore.Qt.LeftButton,QtCore.Qt.LeftButton,QtCore.Qt.NoModifier)
    try:
        self.iface.mapCanvas().mousePressEvent(virtualMouseClick)
        print 'working'
    except:
        print 'not working'

After @sigeal's advice, I have given a go at implementing some things - but run into a bit of a problem.  I am supposing that the sub-classed tool needs its x,y coordinates as pixel rather than map coordinate origin.  Yet use of QgsMapToPixel returns '-inf' (negative infinity) values.

What might be the syntax error I am committing, and am I off target in actions on @sigeal's advice?

#-------From main class----------------------
    ...
    for feature in layerRef.getFeatures():
        if feature.id() != int(featureID):
            pass
        elif feature.id() == int(featureID):
            point = feature.geometry().centroid().asPoint()
            x,y = point[0],point[1]
            deviceX,deviceY = QgsMapToPixel().transform(float(x),float(y))[0],QgsMapToPixel().transform(float(x),float(y))[1]
            print deviceX,deviceY
    self.iface.setActiveLayer(layerRef)
    self.iface.actionIdentify().trigger()
    try:
        featureIdentifyTool(self.iface,layerRef,deviceX,deviceY)
        print 'working'
    except:
        print 'not working'

    ...

class featureIdentifyTool(QgsMapToolIdentifyFeature):
   geomIdentified = pyqtSignal(object, object)
   def __init__(self, iface, layer, x, y):
       QgsMapToolIdentify.__init__(self, iface.mapCanvas())
       self.canvasReleaseEvent(layer,x,y)
   def canvasReleaseEvent(self,layer,x,y):
       results = self.identify(x, y, self.ActiveLayer, [layer], self.VectorLayer)
       print results
       if len(results) > 0:
          self.geomIdentified.emit(results[0].mLayer, results[0].mFeature)


texts: Convert multiple polylines to polygons in QGIS

texts:
I have a vector layer with several polyline segments that represent borders of different areas. Each segment is connected to at least two other segments (at its start and its end). Some segments are connected to more then two other segments.

I need to convert those multiple polylines to multiple polygons. (See fig.)

When I use polygonize tool in vector->geometry menu it treats each segment as it's own polygon and makes garbage. (See fig 2)

Merging all segments into a single one doesn't help. Clipping large polygon with lines also doesn't help in this case since line segments doesn't cross the whole large polygon.

This problem seems obvious but I'm stuck. I can do it in GRASS with ease, but I have no idea how to do it in QGIS. I don't mind using Python scripting for this.


texts: Trouble getting metacarta to work with QGIS

texts:
I have a basic issue with getting QGIS to accept the metacarta WMS service and wonder if anyone can guide me.

I add the "basic layers" to the map, but still can't see the tiles there when I go to the map. I even do the right click -> zoom to layer extent, but there is no visual data in QGIS. Only complete white.

I must be missing something basic so I appreciate any help or advice.


texts: Attribute table : How to select more then one feature?

texts:
I have a vector layer: ROADS that is composed of different types of roads: primary, secondary, tertiary, trunk, etc.

When I right click on the layer ROADS and I go to Attribute table -> Select by expression -> the column where I want to do the selection: highway -> I wrote this:

"highway" ILIKE '%primary%' OR "highway" ILIKE '%secondary%' OR "highway" ILIKE '%tertiary%' OR "highway" ILIKE '%trunk%' NOT "highway" ILIKE '%trunk_link%'

Till NOT "highway" ILIKE '%trunk_link%' the selection is working well. I can see all the roads that starts with primary, secondary, tertiary and trunk. So also trunk_link -> I don't want that. To erase it from the selection I add NOT "highway" ILIKE '%trunk_link%'. It didn't work -> suddenly I had O selected lines .

Do anyone knows the solution to my problem?


texts: is there any way to run a GIS (Esri, Qgis or mapinfo) on a Windows8 RT tablet?

texts: I am looking for a good idea to run a GIS program on a Dell XPS 10 with the Wind.8 RT.Thanks for any


texts: How do you 'incrementally' create vector tiles?

texts:
I am serving my own vector tile layer. In my current setup when a point of interest is added by a user it is stored in PostGis and the relevant zoomtiles are reproduced and written over the previous ones. For high zoomlayers this is fast, however for low zoom (e.g. the world) the processing is demanding as in recreating the vector tiles millions of points of interest are involved.

Is it necessary to 'recalculate' the entire tile? Are there ways to recreate vector tiles in an 'incremental' way which is less processor demanding? This would make sense as the content at lower zooms is >99% the same as the previous vector tile. Does protobuf or other vector tile formats allow such patching? What are existing solutions to my problem?

BTW I am already using a separate tile layer for the points of interest, so it is not like I am recalculating all features. But the points of interest alone can run into millions at world zoom.


texts: How to I change the extents of a Euclidian Distance Raster?

texts: I am calculating the Euclidean distance of a set of points.  The output raster truncates the resulting raster image to the farthest north, south, east, and west point.  I need the distances to extend past those points.  I have created another constant raster that is roughly the size of the area that I need.


texts: Spatial index not speeding up spatialite query containing Multipoint Geometry Data?

texts:
I am querying a sqlite format geodatabase in which there is a label for each geometry shape.

Think of data as :

label = name of a district in your state

geometry: polygon which specify the boundary of that district.

Now, I have coordinates of a person and I want to find out in which state he resides. The problem is that there are around 170 thousand districts and running my query for each input i.e. coordiantes takes around 7 seconds. After reading other articles, I tried spatial index but the results are same. It still takes around the same time. I guess I am not using it properly. I need a bit of advice on how can I improve the query time.

Code involved [Python 3.5]:

voting_district_conn= sqlite3.connect('Voting_Districts.sqlite')
voting_district_conn.enable_load_extension(True)
voting_district_conn.execute("SELECT load_extension('mod_spatialite')")
...
index_cursor = voting_district_conn.execute("SELECT CreateSpatialIndex('voting_districts','GEOMETRY')")
for item in index_cursor:
print(item)  #returning (0,)
...
start_time = time.time()
voting_district_cursor = voting_district_conn.execute("SELECT label from voting_districts where within(GeomFromText('POINT("+lng+" "+lat+")'),GeomFromWKB(GEOMETRY))")
for item in voting_district_cursor:
    print(item)
    print("--- %s seconds ---" % (time.time() - start_time))


texts: Using rasters/features created at 10.0 in ArcGIS Desktop 9.3?

texts:
I am tempted to go back to ArcGIS 9.3 from ArcGIS 10. For what I do, 9.3 is more efficent.

I know I cannot open ArcGIS 10 documents (maps, scenes) and layers in 9.3 unless saved as 9.3 file, but what about features and rasters I created in 10?

Can I add those into any 9.3 document?

Are there any issues that I should expect?


texts:
I am new to gis and am trying to plot a shapefile containing land polygons. (From OSM.) Unfortunately my output image is always empty (i.e. white). Can somebody tell me what I am doing wrong?

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature

fname = r'simplified-land-polygons-complete-3857\simplified_land_polygons.shp'

ax = plt.axes(projection=ccrs.Robinson())
shape_feature = ShapelyFeature(Reader(fname).geometries(),
                                ccrs.PlateCarree(), edgecolor='black')
ax.add_feature(shape_feature)
plt.show()


texts: Seeking examples of Beautiful Maps at county level or closer?

texts:
I was wondering if anyone else had similar examples to Seeking examples of beautiful maps?, but at a county or parcel level.

At those levels, I wonder if people will react the same, or if there will be other facts to consider at that resolution.


texts:
I am trying to clip a layer with a polygon.
Some of my some clippers overlap each other. But it seem to only be the "aggregated" outer line that clip my layer. I want all layer to be clipped by all clippers, i.e. wherever a black line crosses polygon.

The image show an example. As seen on the picture the selected feature is only cut by the outermost border instead of both the elliptical polygons.

To do this I try to use the Clipper with the following parameters