texts:
I'm trying to use a combination of geopandas, pandas and folium to create a polygon map that I can embed  incorporate into a web page.

For some reason, it's not displaying and wonder if anyone can help.

The steps I've taken:

Grabbed a .shp from the UK's OS for Parliamentary boundaries.
I've then used geopandas to change the projection to epsg=4326 and then exported as GeoJSON which takes the following format:

{ "type": "Feature", "properties": { "PCON13CD": "E14000532", "PCON13CDO": "A03", "PCON13NM": "Altrincham and Sale West" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -2.313999519326579, 53.357408280545918 ], [ -2.313941776174758, 53.358341455420039 ], [ -2.31519699483377, 53.359035664493433 ], [ -2.317953152796459, 53.359102954309151 ], [ -2.319855973429864, 53.358581917200119 ],... ] ] ] } },...

Then what I'd like to do is mesh this with a dataframe of constituencies in the following format, dty:

constituency        count
0   Burton          667
1   Cannock Chase   595
2   Cheltenham      22
3   Cheshire East   2
4   Congleton       1
5   Derbyshire Dales    1
6   East Staffordshire  4

import folium
mapf = folium.Map(width=700, height=370, tiles = "Stamen Toner", zoom_start=8, location=     ["53.0219392","-2.1597434"])

mapf.geo_json(geo_path="geo_json_shape2.json",
              data_out="data.json",
              data=dty,
              columns=["constituency","count"],
              key_on="feature.properties.PCON13NM.geometry.type.Polygon",
              fill_color='PuRd',
              fill_opacity=0.7,
              line_opacity=0.2,
              reset="True")

The output from mapf looks like:

mapf.json_data

{'../../Crime_data/staffs_data92.json': [{'Burton': 667,
   'Cannock Chase': 595,
   'Cheltenham': 22,
   'Cheshire East': 2,
   'Congleton': 1,
   'Derbyshire Dales': 1,
   'East Staffordshire': 4,
   'Lichfield': 438,
   'Newcastle-under-Lyme': 543,
   'North Warwickshire': 1,
   'Shropshire': 17,
   'South Staffordshire': 358,
   'Stafford': 623,
   'Staffordshire Moorlands': 359,
   'Stoke-on-Trent Central': 1053,
   'Stoke-on-Trent North': 921,
   'Stoke-on-Trent South': 766,
   'Stone': 270,
   'Tamworth': 600,
   'Walsall': 1}]}

Although the mapf.create_map() function successfully creates a map, the polygons don't render.

Can anyone suggest any debugging steps?

I'm never clear on how to add the full data files if anyone needs them, so please let me know.

Here's the input json file: https://www.dropbox.com/s/9hbfwb95864uiru/geo_json_shape2.json?dl=0

Here's the output data.json:

https://www.dropbox.com/s/x9wh2bodr8feqo9/staffs_data89.json?dl=0


texts: Creating shape around Europe especially France using QGIS?

texts:
I want to create shape around Europe especially for France. But I'm not really sure how to do it as I am new in QGIS. Can I get any good suggestion for this task?

I'm attaching a picture which shows what I actually want. As it is presented yellow line shows the shape which we made using QGIS and I want to do the same thing around France. Any suggestions so far? Is there any external package do I need to use? Moreover I want to do this work automatically rather than manually.

I'm using QGIS 2.18.2 version.


texts: Splitting zAware polyline by points using ArcObjects with C#?

texts:
I'm trying to split a polyline by the boundaries of a polygon using IPolyCurve2.SplitAtPoints using C#.
The polyline I'm trying to split is zAware and holds elevation values at a regular interval.
After splitting the polyline I would obviously like to get the different parts and save them as new features in my feature class, then deleting the original one.
The problem is that when I try to set feature.shape to the geometry in my geometry collection made up of all the polyline parts, I get the error message that the geometry has no Z values.

Does anybody know how to solve this issue? I've looked around in the forum but all I found was code that doesn't deal with polylines that hold Z values.

Here is the code I'm using, passing in the original polyline, the point collection that holds all the intersection points with the polygon and the original feature class I create new features in.

public void SplitPolylineFeature(IFeature pPolylineFeature, IPointCollection pSplitPointCollection, IFeatureClass pFeatureClass)
    {
        //split the feature, each split makes a new part
        IEnumVertex pEnumVertex = pSplitPointCollection.EnumVertices;
        IPolycurve2 pPolyCurve = pPolylineFeature.Shape as IPolycurve2;
        IEnumSplitPoint pEnumSplitPoint = pPolyCurve.SplitAtPoints(pEnumVertex, true, true, -1);
        object Missing = Type.Missing;

        if(pEnumSplitPoint.SplitHappened)
        {
            //new geocoll for polycurve
            IGeometryCollection pGeometryCollection = pPolyCurve as IGeometryCollection;

            //loop through the parts of the split polyline
            for(int intPartCount = 0; intPartCount < pGeometryCollection.GeometryCount; intPartCount++)
            {
                IGeometryCollection pLineGeoColl = new PolylineClass();
                IGeometry pGeometry = pGeometryCollection.get_Geometry(intPartCount);
                //IZAware zAware1 = (IZAware)pGeometry;
                //zAware1.ZAware = true;

                pLineGeoColl.AddGeometry(pGeometry, ref Missing, ref Missing);

                IFeature pFeature = pFeatureClass.CreateFeature();

                pFeature.Shape = pLineGeoColl as IGeometry;  //code crashes here: Geometry has no Z values
                pFeature.Store();

            }
        }
        pPolylineFeature.Delete();
    }


texts: Why does the added raster layer have strange edges?

texts:
I'm using Geoserver in conjunction with the google maps API to add a raster layer onto a map.  It loads fine, however there's some strange problem with opacity at the edge of the layer where you can see a boundary.  Have you any idea how to fix this?

(The original is at obstest.heliohost.org/map2.html.)

I'm totally new to gis and to Geoserver so sorry if this is an obvious question!

Thanks.


texts:
Instead, it is only clipping a rectangle.
I have not used in QGIS in a while so maybe I'm missing a step?
Create shapefile (polygon in this case) > raster > extraction > clip raster by extent > clipping extent (polygon)
Output:

the light grey is the original image, and the slightly darker image is the ''clipped'' image.
Note how it does not match the polygon.
I am using QGIS 3.16.3 on Ubuntu 20.04


texts: Extracting x,y,z values from a postgis linestring to python using ogr

texts:
I am trying to load a postgis linestring into python (using ST_DumpPoints) and extracting x,y,z coordinates, which is working but extremely slow..and if i mean slow, i really mean it (yes, it takes about 5 seconds per row, so imagine how long this is for around 5000 features).
Funnily enough, if i read the same from a points layer (it is originally taken from a gpx-file, so it has track points and a track), the same amount of points is done in a jiffy. I have no idea where the performance loss is coming from and how to resolve it. I rely on this linestring though to calculate a comparing sample. Is it known that linestrings are not be meant for such kind of processing?

The code is as follows:

from osgeo import ogr
import numpy
source = ogr.Open("PG:host=localhost dbname=tracks user=anush")
if source is None:
    print "connection was unsuccessful"
else:
    gFC = source.ExecuteSQL("SELECT st_numpoints(wkb_geometry) from w_130401_aster_comp.w_130401_aster_comp;")
    gFCi =(gFC.GetFeature(0)).GetField(0)
    tP = source.ExecuteSQL("SELECT (ST_DumpPoints(wkb_geometry)).path as path, ST_X(st_astext((ST_DumpPoints(wkb_geometry)).geom)) as X, ST_Y(st_astext((ST_DumpPoints(wkb_geometry)).geom)) as Y, ST_Z(st_astext((ST_DumpPoints(wkb_geometry)).geom)) as Z from w_130401_aster_comp.w_130401_aster_comp;")

    ele_comp = []
    for i in range(1,gFCi):
        f = tP.GetFeature(i)
        z = f.GetFieldAsDouble(3)
        ele_comp.append(z)

If i use the track_points alternative, the results are collectable pretty much right away, i.e.

from osgeo import ogr

source = ogr.Open("PG:host=localhost dbname=tracks user=postgres")
if source is None:
    print "connection was unsuccessful"
else:
    data = source.ExecuteSQL("SELECT ogc_fid, the_geom, ele from w_130401.track_points")

    ele = {}
    for i in range(1,data.GetFeatureCount()):
        f = data.GetFeature(i)
        coords = f.GetGeometryRef()
        x = coords.GetX(); y = coords.GetY()

Another really strange thing is that for the linestring i cannot use an ogr.GetFeatureCount() which returns 0, so i have to pre-run a st_numpoints to initiate the for loop (which returns the number of vertices as integer).

Does anyone know what is going on here?


texts: OpenLayers dynamic rules defintion & clone function filter

texts:
I'm developing on a OpenLayers application which creates dynamic rules for a vector layer using the geostats library. I want to create these rules inside a simple loop through the number of classes. I'm pushing every rule to a array and add this array via stale.addRules(). In my opinion every rule has to be cloned before adding it to the array. The problem is, I am using a OpenLayers.Filter.Function in which the clone() function is not implemented yet (OpenLayers 2.12-rc4).

How can I solve this?

At the moment, only the last rule gets visualized. It seems like the other rules get lost, due to the reference value issue...

My code looks like this:

for (var i = 0; i < numClasses; i++) {
    filter_x = new OpenLayers.Filter.Function({
            evaluate: function(attributes) {
                [...]
            }
        })
    var rule_x = new OpenLayers.Rule({
        filter: filter_x,
        symbolizer: { fillColor: colors[i],
                    fillOpacity: 0.5, strokeColor: "white"}
    });
    //var clone = filter_x.clone(); //returns null! clone is not implemented
    rules.push(rule_x)

}
style.addRules(rules);


texts: Do different ellipsoids have different location of centre?

texts: I wanted to understand how latitude and longitude of a location change by changing the reference ellipsoid. If the location of centre in 2 different ellipsoids is same, and they differ only in the length of minor/major axis, then I don't see how latitude /longitude of a place will differ if any of those 2 ellipsoid is taken as reference.


texts: OpenLayers3 get feature request for more than one point in the same place

texts:
I created a map that brings up a popup overlay on click, the JSFiddle is here. The only issue I have now is if I have two points in the same location it.

As you can see with the Birmingham point, when you click on the point it only brings up the value of the point written in last (on top), ideally I would like it to show the values of any points that are in this location?


texts:
I'm having a weird problem with grids in QGIS 3.12.2.
I have generated a grid to calculate some zonal statistics, and am using that to display a few different things. All good so far. Nice pretty picture. But then I open up QGIS again to fix something, and it's like the attributes have all been scrambled. so instead of showing this clumpy thing like it should

It looks like this scattered mess. It's not just missing 90% of the data, the cells that it's showing are in the wrong spot.
The attribute tables look like this: Everything out of place, numbers going diagonally. It's also showing a "filter" - which might explain some of the gaps, but it is set to show all features so I'm not sure where those are hiding.

If I look at the numbers in the attribute fields, they don't make any sense - it's a generated grid so it should have sensible IDs and left/right/top/bottom info, but that's all been scrambled with the statistic values
There's nothing I've actually done or changed in between it displaying correctly and getting all messed up, and I'd really prefer not to build them again from scratch.
I've tried rebooting QGIS, and reloading the shapefiles, but no dice. It has happened in the past and recovered as spontaneously as it happened, but I'd love to know how to fix it myself.


texts: Highlight polygon on mouse over in Leaflet?

texts:
I am doing a map on Leaflet and I would like to know how to put my polygons highlighted visually in some way when they are covered with a mouse a bit like that: http://www.gistechsolutions.com/leaflet/DEMO/basic/basic_Poly.html

I know you have to use the function highlightFeature but I can't use it.

$(document).ready(function () {
  map = L.map("map").setView([46.227638, 2.213749], 6);
  lControl = L.control.layers().addTo(map);

  fondPlan = L.tileLayer(
    "https://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/{z}/{y}/{x}",
    {
      maxZoom: 17,
    }
  ).addTo(map);
  lControl.addBaseLayer(fondPlan, "Carte topo");

  fondSAT = L.tileLayer("http://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}", {
    maxZoom: 17,
  }).addTo(map);
  lControl.addBaseLayer(fondSAT, "Carte Satellite");

  $.getJSON(
    "geojson/Muron.geojson",
    {},

    function (data) {
      vmuron = L.geoJSON(data, {
        style: function () {
          return {
            color: "#FF0000",
            weight: 2,
            opacity: 0.7,
          };
        },

        onEachFeature: function (feature, layer) {
          contenu =
            "ID: " +
            feature.properties.prefixe +
            feature.properties.section +
            feature.properties.numero +
            "
";
          contenu +=
            "Surface: " + feature.properties.contenance + " m²
";
          layer.bindTooltip(contenu, {
            permanent: false,
            sticky: false,
            opacity: 1,
            direction: "top",
          });
        },
      });

      vmuron.addTo(map);
      lControl.addOverlay(vmuron, "Muron");
    }
  );
});


texts: ArcGIS 10.1 crashes when a button in an addin invokes GDAL DLL through SWIG CSharp wrapper

texts:
I am developing an addin with a button that tries to do coordinate transformation using GDAL OSR. When I create a new CoordinateTransformation object ArcMap crashes and close. If a debug against a debug osr_csharp.dll I can see that CoordinateTransformation object is throwing an exception with this message: Unable to load PROJ.4 library (proj.dll), creation of\nOGRCoordinateTransformation failed. I have checked the addin file unziping it and I can see it contains the proj.dll and all the dependencies that osr_csharp.dll should have.

The rare thing is that the same code executed inside an EXE project works without problem. Where is ArcGIS deploying addins? I cannot see the folder where ArcGIS is supposed to be uncompressing the addin file?

Any ideas on this type of fail?

Thanks in advance


texts: Remove border lines of SRTM tiles rendered as hillshade

texts:
I have downloaded several SRTM DEM raster tiles from https://dwtkns.com/srtm30m/, and created hillshade versions from the tiles. However, upon closer inspection, the hillshade tiles now have border lines that I do not know how to hide or remove (like in the image below, which shows the border lines of the 4 tiles):

I tried checking the Styles options but cannot find a way to hide the lines/borders. When trying to render the tiles as singleband pseudocolor, no lines are shown, so I am suspecting the lines are only generated when doing a hillshade render:

How can I hide/remove these hillshade tile border lines?


texts: Adding drop-down menu to leaflet map to update map/tiles to new location?

texts:
What I want to do is create an interactive choropleth map that has a dropdown menu of various countries so that the map updates when you click on that country.

As an example, here's an interactive choropleth of Ethiopia that I have, along with the code in the index.html file:

http://laurenoldham1202.github.io/magic-ethiopia/
https://github.com/laurenoldham1202/magic-ethiopia

I want to recreate this map so that it has a drop-down menu of about 2 dozen other countries with the same content. The map will act the same way as it does with just Ethiopia, but the user will have options to choose which country they want to view. When they click on that country, the map will update to show that country.

Is there a way to do this inside of Leaflet?


texts: Cesium 3D - Determining the Map Scale of the viewed Globe

texts:
I am working withe Cesium 3D at the moment and want to determine what the equivalent map scale the globe is being viewed at.

I am able to get the viewing altitude in Meters that he Camera is at, but according to this discussion, I would need the eye distance from the surface as well, and I do not understand this.

My ultimate aim is to determine a map scale, convert that to an Openlayers 2 Zoom Level so I can sync the two views (approximately)


texts: Dymaxion/Butterfly alternative projection for QGIS

texts:
I've been reading the Dymaxion and Butterfly projections are not open source, so can not be used with QGIS.

Are there alternatives to these available in QGIS?

I'm looking for a single projection so I can export relatively-accurate SVG shapes of countries around the world when I need them.

It's important these countries are relatively-close to the way they would appear on a globe.


texts: Adding WMS layer (basemap) to ArcGIS Javascript API?

texts:
I am totally new to javascript and using the ArcGIS API for JavaScript 3.19.

I am writing a script to add a map to our gazetteer search web page where the address searched / found can be displayed as a pin on a map.
I managed so far to do it pointing to our Basemap cache on our internal ArcGIS rest services but I have no idea how to replace this cache with a WMS layer that we use as a service which has a much better coverage (whole of UK) - i.e. external url with a username and password in it (sorry - I can't share the url as it is a paid service)

I have looked at the resources and samples on ArcGIS API for JavaScript site but I'm still confused as it seems to refer to layers when our service is a composite of various layers.


texts:
I am curious as to the location of the symbols group files. I have created numerous custom symbols, as well as symbol groups.

I would like to share these with the people I work with.


texts: Getting Python error when trying to open QGIS processing history - Mac QGIS 3.16

texts:
I am using:
QGIS 3.16,
Python 3.9,
GDAL 3.2,
MacOS 10.15.
I've tried a few different QGIS version but keep getting this error when I would like to see the geoprocessor history log by clicking Processing tab -> History:

WARNING    Traceback (most recent call last):
              File "/Applications/QGIS3.16.app/Contents/MacOS/../Resources/python/plugins/processing/ProcessingPlugin.py", line 352, in openHistory
              dlg = HistoryDialog()
              File "/Applications/QGIS3.16.app/Contents/MacOS/../Resources/python/plugins/processing/gui/HistoryDialog.py", line 77, in __init__
              self.fillTree()
              File "/Applications/QGIS3.16.app/Contents/MacOS/../Resources/python/plugins/processing/gui/HistoryDialog.py", line 107, in fillTree
              entries = ProcessingLog.getLogEntries()
              File "/Applications/QGIS3.16.app/Contents/MacOS/../Resources/python/plugins/processing/core/ProcessingLog.py", line 68, in getLogEntries
              lines = f.readlines()
              File "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/encodings/ascii.py", line 26, in decode
              return codecs.ascii_decode(input, self.errors)[0]
             UnicodeDecodeError: 'ascii' codec can't decode byte 0xc3 in position 1268: ordinal not in range(128)

Upon startup of QGIS I get this general information:
 INFO    Programstatus:
             QGIS_PREFIX_PATH miljövariabel:
             Prefix: /Applications/QGIS3.16.app/Contents/MacOS
             Plugin Sökväg: /Applications/QGIS3.16.app/Contents/MacOS/../PlugIns/qgis
             Paketdata Sökväg: /Applications/QGIS3.16.app/Contents/MacOS/../Resources
             Aktivt Temanummer: default
             Aktivt Tema Sökväg: /Applications/QGIS3.16.app/Contents/MacOS/../Resources/resources/themes/default/icons/
             Standardsökväg för Tema: :/images/themes/default/
             SVG Sökväg: /Applications/QGIS3.app/Contents/MacOS/../Resources/svg/
             ' /Users/agh/Library/Application Support/QGIS/QGIS3/profiles/default/svg/
             ' /Applications/QGIS3.18.app/Contents/MacOS/../Resources/svg/
             ' /Applications/QGIS3.16.app/Contents/MacOS/../Resources/svg/
             Användar DB Sökväg: /Applications/QGIS3.16.app/Contents/MacOS/../Resources/resources/qgis.db
             Autentiserings DB Sökväg: /Users/agh/Library/Application Support/QGIS/QGIS3/profiles/default/qgis-auth.db

In the QGIS -> Preferences -> System, there are no old environment path set.
I do find the actual log-file in:
/Users/agh/Library/Application Support/QGIS/QGIS3/profiles/default/processing/processing.log

Can't find the .py files mentioned though, eg HistoryDialog.py .
Any idea how I can see the history log from the actual QGIS application?


texts: Sometimes when the device is pulling coordinates with the GeoPositionAccuracy set to High and the MovementThreshold set to 0, the GeoCoordinate that is captured is further away and not inline with the other coordinates. How can I prevent that or remove that GeoCoordinate from my collection of GeoCoordinates? You'll notice that in the imge below right above the second light grey 'e' and to the right on the route 11 symbol, you'll see a little jog out. The same thing under the first little grey 'e'.


texts: Trying to match delimited text layer to shapefile

texts:
I'm pretty new to QGIS and having a bit of trouble getting a text layer to match my shapefile. I'm trying to plot Canadian cities on a shapefile I have of Canada provinces.

The provinces shapefile is from: http://www5.statcan.gc.ca/access_acces/alternative_alternatif.action?l=eng&dispext=zip&teng=gpr_000b11a_e.zip&k=%20%20%20%2040968&loc=http://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/files-fichiers/gpr_000b11a_e.zip

and the city data I got from: geocoder.ca under there free data downloads.

I tried reading a few articles on this and enabled on the fly CRS transformation. Looking at the provinces shapefile for the metadata it says:

In layer spatial reference system units : xMin,yMin -141.018,41.6814 : xMax,yMax -52.6194,83.1355
Layer Spatial Reference System:
+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

When I plotted the city file I choose to project it using WGS 84.

In layer spatial reference system units : xMin,yMin 0,0 : xMax,yMax 140.87,76.4935
Layer Spatial Reference System:
+proj=longlat +datum=WGS84 +no_defs

and unfortunately my city data doesn't look plotted right. It looked mirrored and there is a point that is way off (guessing the data is wrong). Not sure why it would be mirrored though.

Suggestions?

Since I don't have the rep yet to post images I uploaded an image of what my output looks to skydrive: http://sdrv.ms/YV2Whl


texts: Is anyone indexing ArcGIS Server services on the web?

texts:
In looking at the REST API, it seems like it wouldn't be that hard to write a bot that crawls the web and indexes services (Map Services, GP Services etc.)

Is anyone already doing this?

Update1:
I can find lots of mapservices with simple google search.

I rather doubt google gives a very high ranking to mapservice pages though, so I'd bet there are a lot that are missed.  Anyway, if someone were to crawl all these mapservice pages and assemble an index and expose that index on a site.  The site could allow me to search for  mapservices that are in a particular region and contain certain keywords and add the mapservices to my map.  AFAIK the arcgis.com map gallery isn't built by a web crawler.

Update2
With a funcion named MapServiceHarvester, it sure sounds like Richie Carmichael might be cooking something up.


texts:
I'm trying to cut donut polygons ready to go into a system that can't display them. Basically the exact opposite of this post: How to convert fake donut holes in polygons to real ones?

I want to cut at any point so that the polygon no longer contains an inner ring(s). (Ideally the shortest path but I don't really care.)

Is this possible in QGIS?


texts: Printing vector layer with labels not in cp1252 encoding via GeoServer MapFish printing module?

texts: I'm printing vector layer via GeoServer Printing Module (MapFish Print 2) of GeoServer 2.15. MapFish only supports fonts COURIER, HELVETICA, TIMES_ROMAN in cp1252 encoding for labels. It's very limiting restriction because cp1252 have very little symbols, but let not discuss such a decision. Instead, I want a solution for using symbols outside that charset in features labels, in particulaly cyrillic. I hasn't even idea how to do that.


texts: How to show two different labels of a particular vector layer with the help of QGIS?

texts:
I am working under a government project in Kolkata with QGIS.

I want to show two different labels of a particular vector layer. just like we can show two different labels in Arc GIS with the help of Append option.

Suppose I want to show  Male population Numbers and Female  population numbers on a District vector layer, so how it can be possible? please tell me the option of this Software (QGIS).

Actually i am able to show one label on a layer but how to show two labels?


texts:
I'm trying to create an OpenLayers online map, using QGis/PostgreSQL(PostGIS) to collect data sets and style the map, and then deploying it online.

I previously deployed using PostGIS, GeoServer and GeoWebCache, and learned to use SLD to stlye the map in GeoServer.

Due to the cost and hassle of using a dedicated server to run PostGIS, GeoServer and GeoWebCache (all on Linux), I moved to using map tiles.  This way I could deploy the online map to a standard (cheap) hosting provider.

I started to use the QGis QTiles plugin, which is still considered experimental.  It appears to generate each tile based on features that fall within the bounds of the tile.  This causes large labels that spill out of one tile, but don't necessarily appear in the next tile.

To find an example of the problem, go to: http://thenipissingroad.com/map/?EPSG=4326&coords=-79.6413,45.66620&zoom=14

The town name is "Magnetawan", and is only appearing on one tile cut-off.  If you look around the map at different zoom levels, you'll see many other labels doing this.

Are there other ways to generate tiles than this plugin, or configure the plugin to solve this problem?

I recall using GeoWebCache, and there was a setting (I think) for it to take into acocunt features in neighbouring tiles.  There was a greater computational cost to generate the tiles, but I think it worked better.


texts:
I am trying to write a code that creates the PRACH per sector per band for a site.

My PRACH data contains: "Site name, lat, long, beamwidth, azimuth, and distance"

I got to the point of creating points, but I can't somehow create the Pie-shaped sector that is represented with the variables that my data contains.

Also, if it were to be created, how can I make it loop till it creates 5000 sectors for instance?

You can find an example in the picture.


texts: How to create an intensive map of Russia using QGIS or any other software?

texts:
I am trying to create an intensive (heat) map of Russia. I have an Excel sheet of 36 Russian regions and their corresponding values.

I know that there were similar questions:
Creating map with states by attribute value using QGIS?

However, I do not know which software is better to use. For the QGIS, I cannot find an appropriate shapefile with Russian administrative regions.

Any help would be highly appreciated!


texts: GeoPandas geocode error when using Nominatim

texts:
I am following this lesson about Geocoding.

I have entered:

# Import necessary modules
import pandas as pd

# Filepath
fp = r"addresses.txt"

# Read the data
data = pd.read_csv(fp, sep=';')
data.head()

#geocoding tool
from geopandas.tools import geocode

#nominatim
geo = geocode(data['addr'], provider='nominatim')
geo.head(2)

but get this error:

  AttributeError: module 'geopy.geocoders' has no attribute 'Yahoo'

What is the problem here?

Also, I have an API key for Nominatim, where would I put that? I defined the key:

key = "#########################"

then tried:

geo = geocode(data['addr'], provider='nominatim', key = key)

and

geo = geocode(data['addr'], provider='nominatim', **kwargs = key)

and they didn't work. I have tried to look for solutions online, but they were unclear about how to use the API key.


texts:
I would like to shift a raster by a specified distance in x and y.

I tried to create a world file, as is suggested here, but neither Raster: Conversion: Translate in QGIS nor  GDAL Conversion Translate in the Processing Toolbox created a world file when saved in the suggested formats.
How can I shift raster layer for making cloud-shadow mask?

I tried from the Python console, as is suggested below, but I get an error: AttributeError: 'NoneType' object has no attribute 'SetGeoTransform'.  I am not familiar with the Python console, so I may be doing something wrong.
How to translate (reposition) a .tif raster layer?

The plugin rasmover, which is suggested in several threads, does not appear to exist anymore?


texts: Printing Shapefile spatial Reference Using ArcPy?

texts:
I am trying to read and print out the coordinate system of a shapefile like this

import arcpy
dataset = r"E:\GIS\Centers.shp"
desc = arcpy.Describe(dataset)

extent = desc.extent
spatialReference = desc.spatialRe

name = desc.name
dataType = desc.dataType
print " File Name => " +name +" /  File Type => " + dataType + " /  Coordinate System => " + spatialReference;

but I am getting this error on printing spatialReference part

  Message   File Name   Line    Position     Traceback
           9        TypeError: coercing to Unicode: need string or buffer, geoprocessing spatial reference object found


texts: Nominatim address locator does't give longitude

texts:
I'm trying to geocode my addresses based on this tutorial:

http://www.digital-geography.com/geocoding-addresses-in-arcgis-the-other-approach/#.VzMq_IR97Rb

But Whenever I use Nominatim geocoder it gives me this error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 10, in <module>
AttributeError: 'NoneType' object has no attribute 'longitude'

The process works fine with Google geocoder. Doesn't Nominatim recognize longitudes or what? When I checked the address from my CSV file in OSM, I see that the address is not a point, its a house so a polygon. Can this be a reason it doesn't recognize just one longitude?


texts: Which point of elevation level from Google earth is referring to?

texts:
About the Google Maps elevation API, which point of the elevation level is referring to? Whether to the actual ground level, or the elevation level that one can detect from the aerial (eg: at rainforest it is the height of the canopy tree)?

The document says it is the "surface" of the earth. But I am unsure whether it is really possible to know what is the surface of the earth if the land is covered with thick forest


texts: Is there a way to keep labels and symbols from displaying inside an extent indicator?

texts: I've been working on a utilities map and have been asked to turn off the labels and symbols within the extent indicator on the larger scale map and show them only in the detailed map. Is there a way to do this without converting labels to annotation and making a copy of and editing each individual shapefile specifically for this map?


texts:
I would like to remove the official repository from several QGIS installations.

Removing the Repository from the Settings works so far with:

from qgis.PyQt.QtCore import QSettings
settings=QSettings()

test=settings.allKeys()
for i in test:
    if u'Qgis/plugin-repos/' in i:
        if u'/url' in i:
            if settings.value(i)==u'https://plugins.qgis.org/plugins/plugins.xml':
                repo_key = i.rsplit('/', 1)[0]
                settings.setValue(i,'')
                settings.remove(repo_key);

But as soon as QGIS is started the next time the repository is added to the settings again.

How can I remove it permanently?


texts:
I'm getting some problems using openlayers with TypeScript, there are some errors in the .on use, can you help me giving some tips?

Argument of type '(event: Event) => void' is not assignable to parameter of type 'EventsListenerFunctionType'.
Type '(event: Event) => void' is not assignable to type '(evt: Event) => boolean'.
 Types of parameters 'event' and 'evt' are incompatible.
  Type 'ol.events.Event' is not assignable to type 'ol.render.Event'.
    Property 'vectorContext' is missing in type 'Event'.

PS: I need to use the event.context but the function() only receives l.events.Event type.

My Component

 this.imagery.on("precompose", (event: ol.render.Event) => {});


texts: Obtaining Geocoding results with county value from ArcGIS world geocoder

texts:
The ArcGIS World Geocoder returns pretty good results for our use, with two exceptions that are deal breakers and I can't yet find answers.

There don't appear to be county values consistently returned for certain types of searches.  Notice this search for a postal code:

http://geocode.arcgis.com/arcgis/rest/services/World/GeocodeServer/find?text=80205&f=pjson&outFields=*

The "Subregion" is where I would expect to see "county" listed.  Similarly, I have not once yet seen "Nbrhd" populated, though that would be very useful.  Without the county value though, ArcGIS isn't a solution for us.

Is there a way to get the State/Province/Region code instead of the full name?  Same question goes for Country.


texts: Raster Caculator Qgis 2.6 v/s Gdal_calc.py by Terminal

texts:
I need calculated monthly mean of Snowcover data by MODIS product, and then modify this script gdal_calc.py in this section

# create alphabetic list for storing input layers
AlphaList=["A","B","C","D","E","F","G","H","I","J","K","L","M",
       "N","O","P","Q","R","S","T","U","V","W","X","Y","Z","a","b","c","d","e"]

and wrote this script:

!/bin/bash
echo -n "./gdal_calc.py " > comandos_py.sh
array=(A B C D E F G H I J K L M  N O P Q R S T U V W X Y Z a b c d e)
cont=0
for file in MCD10A1_*.tif ; do
  echo "map: $file"
  letra=${array[$cont]}
  echo -n "-$letra ${file} " >> comandos_py.sh
  cont=$[cont+1]
done
echo -n "--outfile=MEAN_2013_06.tif --calc=\"(A+B+C+D+E+F+G+H+I+J+K+M+N+O+P+Q+R+S+T+U+V+W+X+Y+Z+a+b+c+d)/30\" --NoDataValue=-9999 --type=Float32" >> comandos_py.sh

both running well..... but the output file is incorrect...

in another way! load all file on QGIS 2.6, open RASTER CALCULATOR insert 30 image SNOWCOVER DAILY sum and divided by 30 and the output file IS CORRECT....

can someone explain to me why this happens?


texts: Paste to Original Coordinates problem in AutoCAD Map 2012

texts:
I have a problem copy-pasting from a drawing to another blank drawing to a same coordinate position. Here's what I have done:

Selected the shape I wanted to copy, a simple rectangle
representing land parcel.
Clicked "copy".
Start a new blank using ctrl+N
Clicked dropdown below "paste" and choosed "paste to original coordinates"
At this point Cad used to be automatically zoom to full extent to the shape I was copied, but today it is suddenly not. So I used zoom command and selected Extent.
The result was a random polylines completely different to what I have copied previously, in a random position/coordinates.

I have reinstalled my AutoCAD but no solutions.


texts: Naming exported features from attribute information?

texts:
I am new to Python and in the process of learning. As a self test I am trying to export all of the municipalities of a single county feature layer as their own feature layers. I am trying to name each feature layer based on the municipality label column but Python returns this error "ERROR 000210: Cannot create output" on the last line which is the CopyFeatures_management tool. How can I achieve my desired result?

import arcpy
import os
arcpy.env.overwriteOutput = True
arcpy.env.workspace = r"C:\gis\python_training\municipality iteration test\mun test.gdb"
mun = "Municipality"
dict = {}
with arcpy.da.SearchCursor(mun,["OBJECTID", "MUN_LABEL"]) as cursor:
    for row in cursor:
        dict[row[0]] = row[1]
del cursor
dict2 = {}
for key, value in dict.items():
    if value not in dict2.values():
        dict2[key] = value
arcpy.MakeFeatureLayer_management(mun,"munlyr")
for key in dict2:
    arcpy.SelectLayerByAttribute_management("munlyr", "NEW_SELECTION","OBJECTID = " + str(key))
    arcpy.CopyFeatures_management("munlyr", str(value))


texts: How to export raster from PostGIS enabled DB as GTiff with gdal_translate

texts:
When I want to export a GTiff from my PostGIS DB with this code:

gdal_translate -of GTiff PG":host='localhost' port:'5432' dbname='landslide' user='postgres' password='pass' schema='public' table='dgm' mode='2'" test.tiff

I get the following error:

Illegal filename. character <'> not allowed.Illegal filename. character <'> not allowed.ERROR 4: `PG:host='localhost' port:'5432' dbname='landslide' user='postgres' password='pass' schema='public' table='dgm' mode='2'' does not exist in the file system, and is not recognised as a supported dataset name.

GDALOpen failed - 4 `PG:host='localhost' port:'5432' dbname='landslide' user='postgres' password='pass' schema='public' table='dgm' mode='2'' does not exist in the file system, and is not recognised as a supported dataset name.

I followed the instructions from gdal.

Is the syntax incorrect? I tried to change the quotes, but no success.


texts:
I am working on a script that automatically calculates something for me. To be sure, I checked the result of the calculation with QGIS. Here are the statistics:

STATISTICS PYTHON-CODE:

MAXIMUM = 0.9896427989006
MEAN = -0.21678630311319
MINIMUM = -0.69955635070801
STDDEV = 0.24564237169808

STATISTICS QGIS:

MAXIMUM = 1
MEAN = -0.4204168536468
MINIMUM = -0.69955635070801
STDDEV = 0.17747956958404

I created my script based on this: https://gist.github.com/jgomezdans/1336359. Here it is:

def filter2 (band2, band5):
    ds2 = gdal.Open(band2)
    data2 = ds2.ReadAsArray()
    ds5 = gdal.Open(band5)
    data5 = ds5.ReadAsArray()
    check = np.logical_and(data2 > 0, data5 > 0)
    filter2 = np.where(check, (data2 - data5) / (data2  + data5), 0)
    return filter2

In QGIS I used the raster calculator:

(band2 - band5) / (band2 +  band5)

Bands 2 and 5 originate from a landsat 7 image and are converted to reflectance values (so no DN).

I think I am doing something wrong with the check and np.where lines but when I just calculate the ratio in the script as:

filter2 = (data2 - data5) / (data2 + data5)

the answer is NaN.

I am not sure if I should worry about the differences in the statistics, but I think I should. And also, I'd like to know why I get NaN when I do not use the check.


texts: Rasterizing ungridded XYZ (LiDAR) files?

texts:
I have a large set of ungridded xyz text files (LiDAR point cloud) and want to  export them to any raster format (tiff, asc, gridded xyz...) for further GIS processing.

The source: https://www.opengeodata.nrw.de/produkte/geobasis/dgm/dgm1/

Example:

 403088.00,5713796.00,71.00
 403088.00,5713798.00,70.98
 403090.00,5713796.00,71.02
 403090.00,5713798.00,71.02
 403116.00,5713752.00,70.57
 403116.00,5713754.00,70.57
 403118.00,5713752.00,70.60

The files are ungridded, unsorted, up to 3GB of size and consists of millions of points each.

I don't want to use LASTools because it's somewhat proprietary. I am looking for a quick Open Source GIS way with QGIS, SAGA, GRASS, GDAL, PDAL.


texts: Linking multiple map displays in Internet Map Server (IMS)?

texts:
Has anyone heard of any GIS Internet Map Server (IMS) capability program that will allow multiple displays of data to be geo linked.

For example, when one pans on one view it will automatically relocates the extent on the corresponding linked viewer.


texts: I have some WFS layers in OpenLayers, containing lines. I try to make a tool to split these lines and add a symbol in the breakpoint. Then save them to a PostGIS database. How can I implement it?


texts: Incorrect pixel values after NDSI calculation using ArcMap

texts:
I've got the following problems with my project. My task is complete the NDSI index calculation in the Antartic region. A problem is some pixel values are much lower than should be.
I'm using ArcMap 10.4.1 verison. My data came from Landsat OLI 8 satellite (https://earthexplorer.usgs.gov/) and the DEM data from REMA database (https://www.pgc.umn.edu/data/rema/)
I've read about Topographic correction, but It didn't work succesfully. I didn't get higher pixel values behind the mountainous area. Luckily I found in QGIS the Topographic correction and I tried all of the topographic correction methods. In the final results I got the exact same NDSI values everywhere.
(The input values was different after the topographic correction)
I also tried this toolbox, but the ArcMap crashed while the last correction was running. That happend 3 times in the same point of the progress.( https://hub.arcgis.com/datasets/fd6ed990984945d5a9827418e67f26df )


texts: Adding points in existing profile graph with labels of its heights using ArcPy?

texts:
I have created profiles at multiple locations with lines.

Now I want to add points to that graph,
how do I do that?

Here is the script:

import arcpy
import os
import matplotlib.pyplot as plt
import numpy

#settings for datasources
linesPath = "lines"
rasterPath = "raster"
pointsPath = "points"

# settings for output profiles PNG's
profileFolder = r"\test_arcpy\profile"
profilePrefix = "ArcPy_"
profilePostfix = "_new.png"
fldName = "profile_na"

# standard NoData mapping
NoDataValue = -9999

# describe raster
inDesc = arcpy.Describe(rasterPath)
rasMeanCellHeight = inDesc.MeanCellHeight
rasMeanCellWidth = inDesc.MeanCellWidth

# search cursor
fldShape = "SHAPE@"
flds = [fldShape, fldName]
with arcpy.da.SearchCursor(linesPath, flds) as curs:
    for row in curs:
        feat = row[0]
        profileName = row[1]

        # read extent of feature and create point for lower left corner
        extent = feat.extent
        print "Processing: {0}".format(profileName)
        pntLL = arcpy.Point(extent.XMin, extent.YMin)

        # determine number of rows and cols to extract from raster for this feature
        width = extent.width
        height = extent.height
        cols = int(width / rasMeanCellWidth)
        rows = int(height / rasMeanCellHeight)

        # create Numpy array for extent of feature
        arrNP = arcpy.RasterToNumPyArray(rasterPath, pntLL, cols, rows, NoDataValue)

        # create empty arrays for distance (d) and altitude (z) and NAP line (zv)
        d = []
        z = []
        zv = []

        # loop through polyline and extract a point each meter
        for dist in range(0, int(feat.getLength("PLANAR"))):
            XYpointGeom = feat.positionAlongLine(dist, False)
            XYpoint = XYpointGeom.firstPoint

            # translate XYpoint to  row, col (needs additional checking)
            c = int((XYpoint.X - pntLL.X) / rasMeanCellWidth)
            r2 = int((XYpoint.Y - pntLL.Y) / rasMeanCellHeight)
            r = rows - r2
            if c >= cols:
                c = cols-1
            if r >= rows:
                r = rows-1

            # extract value from raster and handle NoData
            zVal = arrNP[r,c]
            if not zVal == NoDataValue:
                d.append(dist)
                z.append(arrNP[r,c])
                zv.append(0)

        # define output profile filename and figure size settings
        elevPNG = profileFolder + os.sep + profilePrefix + profileName + profilePostfix
        fig = plt.figure(figsize=(10, 3.5)) # inches

        # plot profile, define styles
        plt.plot(d,z,'r',linewidth=0.75)
        plt.plot(d,z,'ro',alpha=0.3, markersize=3)
        plt.plot(d,zv,'k--',linewidth=0.5)
        plt.xlabel('Distance from start')
        plt.ylabel('Elevation')
        plt.title('Profile {0} using Python matplotlib'.format(profileName))

        # change font size
        plt.rcParams.update({'font.size': 8})

        # save graph as PNG
        fig.savefig(elevPNG, dpi=300)

# delete some objects
del arrNP, XYpointGeom, XYpoint
print "PNGs stored in: {0}".format(profileFolder)


texts: Zoom in for selected lands in OpenLayers

texts:
I am trying to make a selection from a search result. I managed to make the selection on the map using the code below.
Now my problem is that when I select one or more buildings with the same style, it does not automatically zoom in on the respective lands. I would also like to zoom in on the selected land.
My code:
var extent = ol.extent.createEmpty ();

var selected_stile = new ol.style.Style({
    fill: new ol.style.Fill({color: makePattern()}),
    stroke: new ol.style.Stroke({
      color: [102,255,51], width: 3
    }),
  text: new ol.style.Text({
    font: 'bold 11px "Open Sans", "Arial Unicode MS", "sans-serif"',
    placement: 'center',
    fill: new ol.style.Fill({
      color: 'white', width: 23
    })
  })
});

imobile.getSource().forEachFeature(function(feature){
  feature.setId(feature.get('text'));

  var array = id_teren.split(",");
  for (i=0;i<array.length;i++){
    var data = "";
    if($('#etc_tarla').is(":checked")){ data = "T:"+feature.get('tarla')+"  "; }
    if($('#etc_parcela').is(":checked")){ data = data+"P:"+feature.get('parcela')+"  "; }
    if($('#etc_suprafata').is(":checked")){ data = data+"S:"+feature.get('suprafata')+"MP.  "; }
    if($('#etc_proprietari').is(":checked")){ data = data+"P:"+feature.get('proprietari')+"  "; }
    if($('#etc_eterra').is(":checked")){ data = data+"C.F.:"+feature.get('e_terra')+"  "; }

     if(feature.get('text') == array[i]){
        console.log('GASITTTTTTTTTTTTTTTTT');
        var feature = imobile.getSource().getFeatureById(array[i]);

        selected_stile.getText().setText(data);
        feature.setStyle(selected_stile);
        ol.extent.extend(extent, feature.getGeometry().getExtent());

     }
  }
});


texts: Why does an empty, non-spatial table get flagged as Is_Simple=False in an SQL query?

texts:
I've created an empty table in an Oracle geodatabase. It was created in ArcCatalog by using the catalog tree >> New >> Table tool. The table is registered with the geodatabase and is not registered as versioned. The only field it has is the default OBJECTID field.

Why does it return False when I run an Is_Simple check? Shouldn't it be returning True, since it doesn't participate in extended geodatabase functionality?

SELECT SDE.GDB_UTIL.Is_Simple('ENG', 'TEST_TABLE')
FROM DUAL;

FALSE

Environment:

ArcGIS Desktop 10.3.1
SDE 10.3.1
Oracle 12c (12.1.0.2.0)

Source SQL:
http://desktop.arcgis.com/en/arcmap/10.3/manage-data/using-sql-with-gdbs/is-simple.htm

Related: What type of data can be edited using SQL?:
http://desktop.arcgis.com/en/arcmap/10.3/manage-data/using-sql-with-gdbs/what-type-of-data-can-be-edited-using-sql.htm

Update 2016-08-23:

To test whether this had something to do with being registered with the geodatabase, I created a separate table that was not registered with the geodatabase (used the Create Unregistered Table geoprocessing tool). When I run the SQL query it returns NOT REGISTERED. In my question above, I'm getting FALSE, as opposed to NOT REGISTERED, so I don't think that's the problem; the table is definitely registered with the geodatabase.


texts: How to standardize raster output from 0 to 100 using raster algebra?

texts:
I am working on generating an integrated moisture index based on the following equation:

IMI = (hill shade * 0.5) + (curvature * 0.15) + (flow accumulation * 0.35)

Iverson et al. 1997

Another source

The equation requires that hillshade, curvature and flow accumulation are standardized from 0 - 100.  How can I accomplish this standardization from 0-100 using raster algebra?


texts: What is the expected margin of error for points plotted on EPSG:3857 projection?

texts: I am using the projection EPSG:3857 to plot out points on the map. Since the Earth is closer to the shape of an ellipsoid instead of a sphere, I am interested to find out what is the range of error if I were to plot some points using this spherical mercator projection. Using examples of prominent points of interest such as the tip of Mount Everest, Statue of Liberty, etc, would be helpful.


texts: Layer added to QGIS TimeManager disappears from the TimeManager panel

texts:
I am trying to use the TimeManager plugin in QGIS 3.16.1. I have successfully used this plugin in the past, but I'm running into some odd problems now that I can't figure out.
I have a geopackage layer with data that looks like this:

fid
statefp
countyfp
geoid
...
date
count

1
36
001
36001
...
2020-12-10
1

I'm adding this layer to TimeManager, selecting the date column for the "Start time", leaving everything else at the defaults. After clicking the "Add" button, I change the "Time Format" for for the layer from "WILL BE INFERRED" to "%Y-%m-%d" to match the format of dates in the layer.
If I then click "OK" to close the settings window, then re-open it, the layer I just added is gone.
There are no errors, either displayed on screen or in the log.
It seems unhappy with the date field in particular:

If I select another field (e.g. fid) as the date field, the layer is added persistently.
If I then change the name of the date field back to date, it disappears.
Date-based filters applied directly to the layer work just fine (i.e., I can set the filter to "date" > '2020-12-01' and it filters as expected).

What's going on here?


texts: How to input python code for custom field edit method in QGIS 2.16? -or - Where is "init function"?

texts:
I am trying to implement the python code suggested as an answer to this question: QGIS Custom Form: “Other” category in combo box replaced by line edit [closed] I have the same goal as that question: to have an edit widget for data entry that starts as a dropdown menu (using value map) then change to a text line edit if the user selects "Other" from the menu. I am experienced with QGIS but new to python coding.

I put the code in a text editor and saved as FillForm.py. Then it seems I should in the layer properties type into the Init Function line edit: FillForm.formOpen. I went to the layer properties and couldn't find the Init function line edit. What should I do with this code to get QGIS to implement it?

My layer general properties look like this:

There isn't a text line edit called "Init function," so I don't know where to type FillForm.formOpen.

Instead I tried inputting the code as an action.

But when I select "Other" from the dropdown menu for the field "SiteName" it doesn't switch to a text edit line like I wanted it to. Now I have the option to run this action from the Identify Results menu

Running the action gives me an error


texts: Speeding up (automating) QGIS georeferencing process?

texts:
I georeference a lot of maps with QGis using a combination of visible street intersections on paper maps and right clicking on those intersections inside of maps.google.com and clicking "what's here?", then copying and pasting the latitude and longitude information into the Georeferencer program by hand.

Does anyone out there have any procedures to speed up their georeferencing efforts?

I've done this so many times I'm a little embarassed as a programmer that I'm doing it in such a pedestrian manner.


texts: What are the differences between WMS, WFS, WCS, WPS?

texts:
By looking at

I do not really understand the differences between WMS, WFS, WCS, and WPS.

Can somebody please explain the differences between them?


texts: Interpreting results of QGIS intepolation and Heat map from point shapefile

texts:
I have taken a step interpolation and heat map of the data point shapefile (building Centroid) with QGIS ver-2.8.1, as shown below.

My question  are:

What is the difference of the two methods?
How to explain the results?(Interpolation and heat Map  of points)


texts:
How to show\hide symbols according scales,

I couldn’t figure out if there is a way to hide symbols depending on scales the same way we do with labels.

For example, in the screenshots below, the “community layer” is categorized according the “classification field. The labels can be turned on\off depending on scales. How the symbols can be turned on\off according scales?


texts: How to have GDAL print layers of GeoPDF AND say which are raster vs vector

texts:
My Objective: I would like to use GDAL to convert a GeoPDF. I want the vector layers as shp files and the raster layers as tif files. I want to do this in a programmatic way.

Edit: In reality, I want to do this with many geospatial PDFs. I'm prototyping the workflow using Python, but it will probably end up being C++. (End Edit)

The Problem: Naturally, the command to convert a vector layer differs from a raster layer. And I don't know (again in a programmatic way) which layers are vector and which are raster.

What I've Tried: First, here is my sample data https://www.terragotech.com/images/pdf/webmap_urbansample.pdf.

gdalinfo webmap_urbansample.pdf -mdd LAYERS

gives the layer names:

...
Metadata (LAYERS):
  LAYER_00_NAME=Layers
  LAYER_01_NAME=Layers.BPS_-_Water_Sources
  LAYER_02_NAME=Layers.BPS_-_Facilities
  LAYER_03_NAME=Layers.BPS_-_Buildings
  LAYER_04_NAME=Layers.Sewerage_Man_Holes
  LAYER_05_NAME=Layers.Sewerage_Pump_Stations
  LAYER_06_NAME=Layers.Water_Points
  LAYER_07_NAME=Layers.Roads
  LAYER_08_NAME=Layers.Sewerage_Jump-Ups
  LAYER_09_NAME=Layers.Sewerage_Lines
  LAYER_10_NAME=Layers.Water_Lines
  LAYER_11_NAME=Layers.Cadastral_Boundaries
  LAYER_12_NAME=Layers.Raster_Images
...

I know to look at the data which are vector and which are raster, but I don't know how to parse this information to know whether to use ogr2ogr or gdal_translate to do the conversion.

Then I thought I could use ogrinfo and just diff all the layers to deduce which ones are raster, but ogrinfo gives me:

...
1: Cadastral Boundaries (Polygon)
2: Water Lines (Line String)
3: Sewerage Lines (Line String)
4: Sewerage Jump-Ups (Line String)
5: Roads
6: Water Points (Point)
7: Sewerage Pump Stations (Point)
8: Sewerage Man Holes (Point)
9: BPS - Buildings (Polygon)
10: BPS - Facilities (Polygon)
11: BPS - Water Sources (Point)

So there's not a one-to-one correspondence with the way these are output.

Does anyone know how to have gdal print the GeoPDF layers and indicate which are raster vs. vector?


texts:
Can someone show me the relevant documentation for relationships between tables?

I'm trying to get the geometry for an entry from the relations table but I'm unable to pick out the right information without this

Details:
schema: pgsnapshot (v0.6)
data-source: http://download.geofabrik.de/europe/great-britain/england/hampshire-latest.osm.pbf
db-load-command: osmosis --read-pbf ./hampshire-latest.osm.pbf --write-pgsql-dump directory=. enableBboxBuilder=yes

As far as I get:

select id, user_id, tags->'name' as name, tags->'admin_level' as level from relations limit 1;

   id   | user_id |   name   | level
--------+---------+----------+-------
 108786 | 2089161 | Waverley | 8

select * from relation_members where member_id = 108786 limit 1;

 relation_id | member_id | member_type | member_role | sequence_id
-------------+-----------+-------------+-------------+-------------
       57582 |    108786 | R           | subarea     |           4

select * from nodes where id = 57582;
 id | version | user_id | tstamp | changeset_id | tags | geom
----+---------+---------+--------+--------------+------+------
(0 rows)

select * from ways where id = 57582;
 id | version | user_id | tstamp | changeset_id | tags | nodes
----+---------+---------+--------+--------------+------+-------
(0 rows)


texts: Download and Browse Sentinel data from AWS

texts:
It's been a while since I downloaded or browsed Sentinel data in AWS from http://sentinel-s2-l1c.s3-website.eu-central-1.amazonaws.com/tiles/

But I can not access this website anymore. Is it closed or is there any AWS link that does the same?


texts: Add simple image layers to MapBox / Leaflet?

texts:
Is there an easy way to add simple image layers (with no, or very limited, GIS data) to Mapbox and/or Leaflet?

Google maps seems to have the ability to add custom layers, simply by attaching the SW and NE corners to points on the image to points on the map

I thought that MapBox or Leaflet would have similar functionality, but I'm struggling so far.

I (presumptuously?) assumed that TileMill might allow me to import a non-GIS labelled image and drag and drop it to the correct position - much as you might align two Photoshop layers. However, I seem to only be able to import GIS datasourced layers, the closest of which to my needs seems to be GeoTIFF.


texts: Splitting OSM data into layers by tags (QGIS modeler)

texts:
I have a problem with using the QGIS modeler to split a downloaded .osm file into layers. What I am trying to achieve is a structure like this:

Traffic (lines)

Main Roads
Roads
Small Roads

Natural (polygons)

Farmland
Forest

Symbols (points)

Church
School

... as well as applying a predefined .qml style to each layer.
What I have tried: I've tried installing the QuickOSM plugin and using its functions in the QGIS modeler. I wanted to split the (offline) OSM file into layers (with attributes for further processing) (see image):

What I get when executing this (simple) command is this:

An error has occured while executing Python code:
TypeError: cannot deepcopy this pattern object
Python version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
QGIS version: 2.12.2-Lyon Lyon, b403532

Am I doing something wrong using the QGIS modeler? To me this operation seems logical: Take the OSM file and split it into points, polygons, multipolygons and lines.
Is there a workaround for this bug?
I have submitted a bug report (https://github.com/3liz/QuickOSM/issues/78), but in the meantime I thought of using a SpatiaLite DB as a solution. Now my question is:

How can I create layers from queries to the SpatiaLite DB inside the QGIS modeler? I know SQL, but which tables should I look at (I need the tags as well as the geometry) and which algorithm should I use? (I do not know what every QGIS algorithm does, this is why I tried the QuickOSM plugin, since it should do exactly this for me)


texts:
Using QGIS 2.6.1 Brighton on an IMac (10.8.5). QGIS 1.8 Lisboa with exactly the same two vector shapefiles and the same CSV file works. I joined the file containing the polygons showing the features on my archaeological site to the csv file. Used the query function and in 1.8 all of the features stay the same color except those that match the criteria. They turn another color. I can then save that layer.

In 2.6.1, the same procedure results in my being able to individually select a single feature and all the rest disappear. Presumably that's a new feature in this version and that's fine. But if I use any of the database material in the csv file as a search criteria, it tells me nobody's home, which is obviously not the case. I've tried on several different headings with the same result.

I've completely uninstalled QGIS and all of the associated files, downloaded the latest and still have the osgeo problem of not loading despite it saying that it references 1.11 as it should.

Any ideas on how to fix this?


texts:
Apologies for cross-posting - I have also asked this on StackOverflow...

I am initially pushing MapBox Marker objects to an array via:

var el = document.createElement('div' + index);
el.className = 'marker';
deviceMarkers.push(new mapboxgl.Marker(el, { offset: [-50 / 2, -50 / 2] }).setLngLat([device.lat, device.lon]).addTo(map));

Elsewhere in the code, I extract the marker via:

var deviceMarker = deviceMarkers[index];

I would like to be able to get 'el' from deviceMarker object, in order to alter the icon size, orientation etc. dynamically. How can I extract 'el' as the original div?


texts: Performing isochrone in PostGIS by pgRouting?

texts:
In pgr_drivingDistance for isochrones the person who answered, missed to say the explanation of the answer about what function used by him to get the output for isochrone. I want to do the same but I got the same result like network table and no difference. with this image, I have attached the code with this.

create table htokdisdri as
    SELECT t1.seq,
           t1.id1 AS Node,
           t1.id2 AS Edge,
           t1.cost,
           t2.geom AS edge_geom
    FROM  PGR_DrivingDistance(
           'SELECT id, source, target, traveltime as cost FROM muensterstreets_t',
          (SELECT id
           FROM muensterlocation
           WHERE "Name" = 'HILTRUP'),
          05, false) t1
    Left join   muensterstreets_t t2 on t1.id2 = t2.id
    order by t1.seq;


texts:
I am trying to find what points (buildings) along a coastline are directly visible from a specific stretch of beach.
I have various digital elevation and surface models available and vectors for the lots/building heights that are my observer points and I am able to run viewsheds to create rasters with their field of view no problem. Another bit of complexity is the points represent future development, so I'm trying to find out whether new development built to the maximum legal height will have an impact on the beach, and there are multiple rows of potential buildings.
However, how do I select the observer points that directly view the beach after I have created the viewshed (I have a polygon vector for the beach) and remove any observer points that don't?
From my understanding all the viewshed does is show how many observers see a certain pixel, it doesn't specify which observers.
Is there any way to glean this information?
The end goal is to create an overlay/polygon layer for building zoning conditions. So if a proposed building is going on a lot that is viewable from the beach it must have special conditions etc.
I may be able to do it with trial and error, running a new viewshed every time I remove an observer point, but there are hundreds of points and it's not really feasible or accurate.
I have also considered running the viewshed from the beach but as I'm working with proposed buildings there becomes an issue where if I create hypothetical buildings on the front lots for the viewshed to see or not, they block sight of any potential buildings behind them. so it's not overly helpful.
I'm starting to think that maybe I need to tackle this from a totally different plugin than viewshed.
I just found this on ArcMap that does exactly what I need, is there anything similar on QGIS?
https://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/observer-points.htm


texts: Split polylines at a specific distance in QGIS / ogr2ogr

texts:
Similar to this question I have a set of polylines that I need to be split by a equal distance. How do I do this with QGIS, gdal/ogr2ogr and other open source tools (on Mac)? I don't have access to ArcGIS, therefore not a duplicate question, I hope.

I figured out to extract the coordinates at a specific distance from start of a line with this in QGIS (1000 m from start):

 x(line_interpolate_point($geometry, 1000))
 y(line_interpolate_point($geometry, 1000))

From there I'm not sure how to continue.

I need this to be able to create points along contour lines where I can extract angle and elevation value (from attribute). See this question and answer where I need this information.


texts: Converting GeoTIFF to USGS DEM or DTED using GDAL

texts:
I have been looking the documentation for gdal_translate and it doesn't list the -of format available options there. So I tried looking at the source and  basically ended up read the SetFromUserInput function source looks for several strings none of which has any reference to dem or dted. For reference here is the link:

https://github.com/OSGeo/gdal/blob/501ca3703ff4d33badbc539ab7121f6263ff1ac3/gdal/ogr/ogrspatialreference.cpp

I see alot of examples like so:

gdal_translate -of KMLSUPEROVERLAY mosaic.vrt out.kmz

but can't find any info on the type of output formats. is there a proper documentation for this?


texts: How does ArcGIS Desktop calculate INSIDE point for Feature to Point?

texts:
Let me back up and say that from what I understand there are two common ways to calculate the centroid of a polygon in ArcGIS Desktop:

Using Calculate Geometry on fields within the attribute table of the feature class.
Using Data Management -> Features -> Feature to Point from the toolbox.

These both give the same result - the geometric centroid of the polygon. However, there is no guarantee that point lies inside the polygon.

The Feature to Point tool has an inside checkbox option, that according to the documentation:

  Uses a location contained by an input
  feature as its output point location.

What I would like to know is how is this point calculated by ArcGIS Desktop and what is its "theoretical" meaning, if that makes sense.


texts:
We have a flat file of 50 Geology maps that have a transparent overlay with faults and drill holes.  I have them scanned and georeferenced as TIF images and they will be used in ArcGIS 10.0, 10.1 and Geosoft Oasis Montaj.

To scan the overlay (which is mostly transparent with black circles and lines depicting the features) I needed to tape it to a white sheet, otherwise the scanner wouldn't recognize that there was anything there.  As a result, the image has a light gray background.  To see the Geology map and the overlay together, the overlay can be displayed as transparent, but that washes out the color of the features too (they turn a lighter grey instead of black).

I don't konw how to work well with raster images, so could someone tell me how I can go about turning light gray background of the overlay to a different value?  I want it to show up as nodata in Arc, but I'm not sure how to grab a range of cell values and change them to a different value.  Thanks.

btw, I have access to Spatial Analyst as well, but I don't know enough about rasters to know where to start.


texts: Area Difference - Attribute Table and Map

texts: How can I setup plot areas to be identical in attribute table and map. In map and in attribute table I'm using the same " $area " designator. What's more interesting, the area of smaller plot is higher than the bigger one..


texts: Converting "Degrees Decimal Minutes" to "Decimal degrees" in Python

texts:
I'm trying to convert the following format to decimal degrees: 45° 17,896' N
Is there any module that I can use? I tried using dms2dec:
from dms2dec.dms_convert import dms2dec
lat = dms2dec("45° 17,896' N")
print(lat)

The output is not correct (position is somewhat wrong), so I assume that dms2dec does not support "Degrees Decimal Minutes" format. Any advice?


texts: Creating random points using QGIS using Random points tool?

texts:
I have a polygon layer that has plot layouts of an area.

I want to create a random point on each polygon. I tried to use the Random points tool (Vector->Research Tool->Random Points) to create the points but its not creating the points directly over the layer. The polygon layer has 1819 polygons. The tool creates the same number of points but at a different location and a different scale.  I could have moved the points layer to the right place but the scales do not match so that method would not work. Can anyone help me out or point me in the right direction.


texts: How to use a conditional logic on Impervious Mapping

texts:
I am trying to use a logic ( rule -based classification ) to extract built-up area, similar to as suggested by Xu (2007).He uses, three spectral indices, SAVI( Soil Adjusted Vegetation Index), NDBI ( Normalized Difference Built-up Index ) and MNDWI ( Modified Normalized Difference Water Index) to help extract built-up area. However, I am not able to know the method he chose a threshold of SAVI=-0.344 in the following statement:( In other words, how am i able to extract the right threshold value for SAVI to avoid spectral confusion?)

Given, Band 1=SAVI , Band 2= NDBI and Band 3=MNDWI

If Band 1 <  -0.344 and Band 2 > Band 3 then 1 Else 0.

  The maximum of built-up land class in band1 (SAVIband) of Fuzhou image
  is -0.343, whereas the minimum of vegetation class in that band is
  - 0.182. Therefore, using -0.344 as a threshold value can help avoid
  the confusion between vegetation and built-up land classes and greatly
  increase the extraction accuracy (see sections below).

Xu, H., 2007. Extraction of Urban Built-up Land Features from Landsat Imagery Using a Thematic oriented Index Combination Technique. Photogrammetric Engineering & Remote Sensing 73 12: 1381–1391.


texts:
Does anyone know, what heuristic function is used in a*?

May i change that heuristic in a* ? and how i change it?

Here some heuristic,
Euclidean distance,
Manhattan distance,
Maximum distance

Thanks :)


texts: Converting map consisting of multiple DXF layers to shape in QGIS 3.10 LTR

texts:
I'd like to develop a workflow for the following issue:
There is a map in DXF format which consists of multiple (924) layers. All together they show a complex city map.
Using the Import/Export function works fine for me. Is there a workflow to import the DXF files and build a single vector layer containing all the DXF layers. A shapefile would be fine for me.


texts:
I am trying to get cloud-free imagery for a specific region during specific dates. I am very new to GEE, but the script below seems to work. However when I look at the imagery that is then printed on the map, all images are black. I've checked and the RGB is correct.

What am I doing wrong?

// A rectangle representing the roi.
var geometry = ee.Geometry.Rectangle(-8.372, 10.017, -2.703, 5.358)

var ic = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
.filterBounds(geometry)
.filterDate('2013-01-01','2013-12-01');

    var c = ic.filterBounds(geometry);

    var withCloudiness = c.map(function(image) {
      var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
      var cloudiness = cloud.reduceRegion({
        reducer: 'mean',
        geometry: geometry,
        scale: 30,
        maxPixels: 1e9
      });
      return image.set(cloudiness);
    });

var filteredCollection = withCloudiness.filter(ee.Filter.lt('cloud', 10));

// Compute the median in each band, in each pixel.
//var median = filteredCollection.median();

print(filteredCollection);

// Select the red, green and blue bands.
//var filteredCollection = median.select('B3', 'B2', 'B1');
var visParams = {bands: ['B4', 'B3', 'B2'], max:0.4};
Map.addLayer(filteredCollection);
Map.setCenter(-6.01, 6.51, 5);


texts: QGIS Server automatically includes port numbers in the GetCapabilities response URLs. How can I change this?

texts:
I have a QGIS server up and running, but my GetCapabilities response always serves all URLs with port number included. Problem is that these URLs only work without port number. How can I make sure the port number is not used in the GetCapabilities request? I am not sure if this is something I should fix in a config file from QGIS Server, or maybe in the web server configuration.

I have a lighttpd web server which is exposed on port 8080. However if I expose this on port 80 (default from lighttpd) my GetCapabilities will not show the port number. Problem is that we are not allowed to expose port 80, so this is not an option.

(Maybe) some relevant information:
Lighttpd.conf:

server.document-root = "/"
server.port = 8080

server.modules += ( "mod_setenv" )
server.modules += ( "mod_fastcgi" )
server.modules += ( "mod_accesslog" )
server.modules += ( "mod_rewrite" )

server.errorlog = "/dev/stderr"
accesslog.filename = "/dev/stderr"

setenv.add-environment = ("QGIS_PROJECT_FILE" => "/geo-map/source.qgs")

fastcgi.server = (
  "/" => (
    "mapserver" => (
      "socket" => "/tmp/mapserver-fastcgi.socket",
      "check-local" => "disable",
      "bin-path" => "/usr/libexec/qgis/qgis_mapserv.fcgi",
      "min-procs" => env.MIN_PROCS,
      "max-procs" => env.MAX_PROCS,
      "max-load-per-proc" => env.MAX_LOAD_PER_PROC,
      "idle-timeout" => env.IDLE_TIMEOUT
    )
  )
)

here you see part of the GetLegendGraphic where 8080 is generated, but the link only works without the port number.


texts: ArcObjects equivalent to ArcMap Flip button?

texts:
In an ArcMap extension, I need to flip multiple lines from a selection.

The ArcMap Flip button only allows you to flip one feature at a time.

What's the ArcObjects equivalent to ArcMap's Flip button?


texts: FME metafile multi-writer could not be opened

texts:
I run a workflow and FME showed me that message :

  2019-08-23 12:46:05|   1.9|  0.0|ERROR |BADNEWS: File C:\Program
  Files\FME\metafile\MULTI_WRITER.fmf could not be opened
  (tabrdr.cpp:161)

It is the very first time I come across that, I don't understand why


texts: applyEdits to FeatureLayer does not add geometry


texts: Solving computing time problem when grid size element increases in flood model?

texts:
I want to create high resolution and high accuracy "area of inundation" map using FLO-2D Model. I have high resolution Digital Elevation Model (DEM) and I must show small topographic details on my maps. To avoid the loss of these details, I need to use small grid size nearly 30 centimeters. If I select small grid size, grid number increases. So, my insufficient computer doesn't solve numerical schema of flood/rainfall modelling.

Do I divide study area some segments?

Or what can I solve this problem?


texts:
I have a shapefile that I have posted to geoserver using the REST api.  I have set the coordinate system to EPSG:3785 as indicated below.

When looking at geoserver this is what I get when looking at the posted layer:

I'm wondering why the declared srs is not being populated here even though it is set on the shapefile and in the list in geoserver.  I've tried several different projections that are in that list that do populate.


texts: Selecting Layer by Attribute using SQL and variable in ArcPy

texts:
I've had a good look around the site and haven't been able able to find any 'select by attribute' questions regarding string variables and how to integrate them into an SQL statement. In my particular case I have used a search cursor to create a variable but cannot find a way to use that variable in the select by attributes tool. Below is an example of one of the many ways I have tried to select the variable.

The variable I've mentioned (point_id) is the string EDRN_044851.

rows = arcpy.SearchCursor("feature_layer11")
for row in rows:
    point_id = row.getValue("DRN_update")

arcpy.MakeFeatureLayer_management("EDRN_NODE", "feature_layer15")
arcpy.SelectLayerByAttribute_management ("feature_layer15", "NEW_SELECTION", '\"LCC_DRN_ID\" = {0}'.format(point_id))

desc1 = arcpy.Describe("feature_layer15")
shapefieldname1 = desc1.ShapeFieldName
rows = arcpy.SearchCursor("feature_layer15")
for row in rows:
    feat = row.getValue(shapefieldname1)
    pnt = feat.getPart()
    pnt_x = pnt.X
    pnt_y = pnt.Y


texts: Cloud-masked MODIS Albedo in Google Earth Engine

texts:
I would like to extract the Albedo values at a certain point. I found an answer in a Stackoverflow post to help my problem and adjusted it. However, I cannot manage to get the cloud-masking to work because I cannot find the correct band to choose for the mask. It keeps showing an error.
And I was wondering if there might be a possibility to include the shape coordinates of the MODIS pixel that is used at the location?
// A function to mask out cloudy pixels.
var maskClouds = function(image) {
  // Select the QA band.
  var QA = image.select('BRDF_Albedo_Band_Mandatory_Quality_shortwave Bitmask')
  // Make a mask to get bit 10, the internal_cloud_algorithm_flag bit.
  var bitMask = 1 << 10;
  // Return an image masking out cloudy areas.
  return image.updateMask(QA.bitwiseAnd(bitMask).eq(0))
}

var POI = ee.Geometry.Point([-49.31582,69.56833]);
var modis = ee.ImageCollection('MODIS/006/MCD43A3');
var modisAlbedo = modis.filterBounds(POI)
                    .filterDate('2019-01-01', '2019-12-31')
                    .map(maskClouds)
                    .select('Albedo_WSA_shortwave');

// Scaling
modisAlbedo = modisAlbedo.map(function(img){
  var date = img.get('system:time_start');
  return img.multiply(0.001).set('system_time_start', date);
});

// Create a function that takes an image, calculates the mean over a
// geometry and returns the value and the corresponding date as a
// feature.
var createAlbedo = function(img){
var date = img.get('system_time_start');
var value = img.reduceRegion(ee.Reducer.first(), POI).get('Albedo_WSA_shortwave');
  var ft = ee.Feature(null, {'system:time_start': date,
                             'date': ee.Date(date).format('Y/M/d'),
                             'value': value});
  return ft;
};

// Apply the function to each image
var Albedo = modisAlbedo.map(createAlbedo);

// Export the time-series as a csv.
Export.table.toDrive({collection: Albedo,
  description: 'MODIStest',
  folder: 'MODIS_Albedo_006_MCD43A3',
  selectors: 'date, value'});

print (Albedo);


texts: What command to use in PostGIS to get same result as ArcGIS?

texts:
I have two shapefiles and i'm trying to do simple spatial analysis. When I want to know which polygon belongs to other polygon in other shapefile I always use select by location and in Spatial selection method i choose contain the source layer feature.

Now i have these two shapefiles in postgis and i'm trying to do the same using SQL Query but i got these results:

Select * from state where ST_Contains(geom, (select geom from region where REGION_N='CAPITAL_REGION'));

=0 results, in ArcMap i got 1 result.

I also tried

  ST_ContainsProperly

=0 results

When i tried

  ST_Intersects

= 3 results, the same in ArcMap. So it looks like intersects works but "ST_Contains" doesn't work for some reason.

Where could be with SQL problem?


texts: Improve speed of PostGIS "nearest neighbor" query?

texts:
One of the features of the Python/Django/PostgreSQL web application I'm building is that a user will be able to see the 100 other users who are nearest to them, sorted by distance.  I'm using PostGIS 2.1 and its "<->" operator for doing the distance ordering and nearest neighbor limits.  However, I'm wondering if I can make my query any faster since this will probably be the most often-run query in my application.

Here's the table:

CREATE TABLE user_account(
    id SERIAL PRIMARY KEY,
    name VARCHAR(64),
    address VARCHAR(64),
    city VARCHAR(24),
    region VARCHAR(18),
    postal_code VARCHAR(10),
    country VARCHAR(2),
    lon DOUBLE PRECISION NOT NULL CHECK(lon > -180 and lon <= 180),
    lat DOUBLE PRECISION NOT NULL CHECK(lat > -90 and lat <= 90),
    location GEOMETRY(POINT, 4326)  -- PostGIS geom field with SRID 4326
);

I've also created an index on the location geometry field:

CREATE INDEX user_account_idx ON user_account USING GIST(location);
VACUUM ANALYZE user_account (location);

Here's the query:

EXPLAIN ANALYZE
SELECT name, city,
       ST_Distance(location::geography, 'SRID=4326;POINT(-118.240584 34.072918)'::geography)/1000 as dist_km
FROM user_account
ORDER BY location <-> 'SRID=4326;POINT(-118.240584 34.072918)' limit 100;

Here's the explanation:

 Limit  (cost=55.55..55.80 rows=100 width=244) (actual time=0.351..0.380 rows=7 loops=1)
   ->  Sort  (cost=55.55..55.92 rows=150 width=244) (actual time=0.347..0.356 rows=7 loops=1)
         Sort Key: ((location <-> '0101000020E610000037FC6EBA658F5DC016A5846055094140'::geometry))
         Sort Method: quicksort  Memory: 25kB
         ->  Seq Scan on user_account  (cost=0.00..50.12 rows=150 width=244) (actual time=0.235..0.313 rows=7 loops=1)
 Total runtime: 0.430 ms
(6 rows)

The reason I'm concerned about speed is that when I ran EXPLAIN ANALYZE on this query with only seven rows in my table, it took a little over .400 milliseconds.  It would appear that having to cast the result to the "geography" type is having a big impact because the time drops to around .185 milliseconds if I omit the cast and get the results in degrees.  But degrees aren't useful when you need to calculate a distance.  Can I make this query any faster while still calculating a distance?  I've read similar questions here and didn't see anything useful.

Thanks for your advice!


texts: Unable to adjust raster symbology in ArcGIS

texts:
I am bringing in raster (tif) data (ocean pH) into ArcGIS (exported from QGIS).
I expected it to come in in a way that I could adjust the colours.  However, it is coming in as RGB colour bands (three bands) and I am unable to query the pH values. Is there a different way to bring in the data so I would be able to adjust the colours and query the values?


texts: Automating import of multiple rasters as layers into tilemill?

texts:
How do I automate the import of raster tiffs into tile mill?

I have well over 100 files.

I have many warped tiff files ready to be imported as individual layers into tile mill so they can be saved out as a single mbtiles file.


texts: Selecting high altitude points for training samples

texts:
I am creating training samples for land cover classes in ArcGIS 10.3, and I need to define mountain areas as land that is over a certain height.

How can I do this?

(I have tried drawing polygons over the areas color coded to be over 600m altitude in my DEM, but I worry that this is irrelevant as the spectral signatures of these mountain areas will be similar to other land cover classes)


texts: Updating a column in MapInfo based on Symbology or colour

texts:
I am transferring a large dataset from MapInfo (V 12.0.2) to QGIS (V2.8.1) and I want to be able to  create a column in the MapInfo attribute table that differentiates features with different symbology (even if it differentiates with a number or letter), that way when I add the data to QGIS I can easily change the symbology to categorized style and re-format them the way they were in MapInfo. I styled my symbols using the style function, I did not use the thematic map tool.

I was able to do this with text data - meaning I could populate the column with the text that had been stored in the file, I am just wondering if this can be done with symbols or polygons?

I have tried the MapInfo-to-QGIS Style generator plugin and it would not work for me.


texts: Delete (don't hide) text annotations in QGIS

texts:
I feel like this question should have a really simple answer, but I can't find it anywhere.

I added an annotation to my map while experimenting with different buttons, and it wasn't at all what I wanted, so I tried to figure out how to delete it. All I could find was the instruction to hit Ctrl + T to toggle visibility off, so I did that. The next day, while again trying to solve a different labeling problem, I hit the text annotation button again, but again, it wasn't at all what I wanted, so I tried to delete. But now I'm stuck because anytime I hit Ctrl + T, either the first annotation I accidentally created is visible, or the second one is.

How can I fully delete these annotations? There must be a way, right? I'm not stuck forever with these unwanted call-out boxes?


texts: How to create a GRID with customized Labels in QGIS

texts:
I would like to know how can I create a GRID layer with customized labels as you can see in the image below:

I have searched this topic here:

Creating indexed vector grid in QGIS?
QGIS Python Plugins Repository
Create a grid with all polygons labelled - Index style

I did not found a clear idea, first of all I activate the Experimental plugins to install "Create Indexed Vector Grid" plugins but it does not appears.
I went to the repository to understand it but when I add the repository it gives this error:

A part from this I used 2 different ways off creating:

1st - Vector > Research Tools > Create GRID...

2nd - MMQGIS > Create > Create GRID Layer...

As you can see the result it the same.

How can I do this?

I´m using QGIS 3.6


texts: Overlay a Raster on an Interactive Map in Python

texts:
I would like to know how I can overlay any raster created from a satellite product over an interactive map. This raster can be a single layer raster like a classification map or a RGB as shown below.

which I would like to overlay with a colourmap like this

I found some solutions for using Jupyter notebook but I would like to have these maps in Spyder.


texts:
For example, I have georeferenced layer loaded and I want lat/lot grid or at least some ruler in QGIS display canvas, to suite me as a helper.
It looks like only reference is cursor tracking in status bar, which is not what I'm after.

Is there such option, or a plugin?


texts: Polygon vertices disappearing/merging when less than 5 ft apart in ArcGIS Desktop?

texts:
I am having trouble with ArcGIS 10.2.2 on a database I received from a client.  Anytime I drawn a new polygon, if the vertices are closer than about 5 feet they will merge into one vertex at the midpoint of the line segment between the two vertices when completing the sketch (F2).  I have looked into possible topology issues but cannot find that topology has been enabled on the database.  The XY Tolerance is set to 0.0032 US Foot.  Snapping is set to 15 pixels and is working.

The sketch line looks fine with vertices at all locations needed.  Press F2 and the final polygon is changed as noted above.  No issues with lines that have vertices closer than 5 feet, just polygons.

Any suggestions?


texts: Zooming to selection in arcpy.mapping and only seeing selected features?

texts:
I selected data by using select layer by attribute, then I need to only select these selection and zoom in to the selection.
my code is this, it can zoom in base on the scale I set, but I see other adjacent polygon that has other Code_No, why?

Layer = arcpy.mapping.ListLayers(mxd, "lyr", df)
arcpy.SelectLayerByAttribute_management(lyr, "NEW_SELECTION", "Code_NO = '%s'"%mycodeNo)
df.extent = lyrLayer1.getSelectedExtent()
df.scale *= 2


texts:
I created a point shapefile by converting a raster into points. The default fields generated were FID, Shape, POINT_ID, and GRID_CODE.

I added two new fields "X_location" and "Y_location". I generated the values for these two newly added fields by calculating geometry. Further, I extracted the corresponding NDVI and LST values for the points in the shapefile.

Now the problem is, when I apply the "Select by attributes" option, I can select values for the fields POINT_ID or GRID_CODE. But when the same is applied for the X_location, Y_location, NDVI or LST field, it doesn't work.

More specifically saying, I applied the following query.

Select by attributes--> Chose method as "Create new selection"--> "POINTID"--> Get Unique values--> "POINTID"= (Some value chosen from the retrieved value list)--> Apply. This query worked. But when the same was applied with a different field like "X_location", "Y_location", NDVI or LST, it selected nothing! Even after selecting a value from the list (after clicking "get unique values"), it showed, "0 out of 2799326 selected". I was only changing the field and getting the new unique values for the same field!

I made sure that the value I wanted to select, was present in the table. Also, in the ArcCatalog, in the file properties, I checked for index and there any index was not present. So I made new ones. Still it did not help.

Where exactly I am doing wrong?


texts: Checking for white space in fields of Excel spreadsheets?

texts:
I just wanted to ask this for future reference, is there a good reliable way to check for white space in excel spreadsheets and such you receive from others?

Three or four times now I've had problems performing joins because of this, I just wanted to know if there is a good consistent way to check for it.


texts: Writing a number if another field is not null in QGIS Table

texts:
I am trying to make a condition in a new column that if there is a value in another column, a number will be filled accordingly
For example
Here's the table

If there is a value in "watt_3" I want to have 3
If there's only in "watt_1" I want to get the value of 1
Iv'e thought of these logical rule because I not really good at coding, it's not correct
CASE

WHEN "watt_3"  NOT is_empty_or_null THEN 3

ELSE

WHEN "watt_2"  NOT is_empty_or_null THEN 2

ELSE 1

END

How can I solve this?


texts: Merge spatial and non-spatial data and create SpatialPolygonsDataFrame in r

texts:
I am trying to merge non-spatial data (data frame in R) and spatial data (SpatialPolygonsDataFrame in R) and eventually create the merged file in SpatialPolygonsDataFrame form.

In specific, I downloaded 'zip areas boundary file' from census web site (https://www.census.gov/geo/maps-data/data/cbf/cbf_zcta.html) and read it as 'SpatialPolygonsDataFrame' file using 'readShapePoly' function in R.

And then I merged the SpatialPolygonsDataFrame with my data in 'data frame' form, but the merged file is 'data frame' form not 'SpatialPolygonsDataFrame'.

Can someone let me know how to create 'SpatialPolygonsDataFrame' after merging spatial data and non-spatial data? I have been spending my entire day just on this but didn't have clues.

The code that I used is as follows:

nation <- readShapePoly("C:/Users/Research/data/nation.shp")
us_urban_zipcode <- as.data.frame(us_urban_index13[c(39,91)])
us_urban_index_nation <-merge(us_urban_zipcode, nation, by = "ZCTA5CE10",
all.x=TRUE)

Nation is SpatialPolygonsDataFrame and us_urban_zipcode is dataframe, and merging them results in dataframe, not SpatialPolygonsDataFrame, which I need for further analysis.

My data non-spatial data looks like this:

---------------
zipcode | row
 10003  |  1
 10002  |  2
 10003  |  3
 10004  |  4
 10002  |  5
  ...   | ..
-----------------

And my spatial data looks like this:

----------------------
zipcode | AFFGEOID10
 10001  |  477175
 10002  |  2118827
 10003  |  78588518
 10004  |  9804480
 10005  |  4845242
  ...   |    ..
-----------------------

So, basically my non-spatial data is bigger than spatial data in terms of observations. The zipcode of spatial data is all unique (only one zipcode in each observation), but there are redundant zipcodes for non-spatial data.
But I need to keep all observation in non-spatial data for further analysis.
This is why I used 'all.x =T' or 'all.y = T' in the merge function.


texts: Line along points with no continous ID and creation date/time

texts:
I have a huge points file. The points are along a route network. I would like to build lines by going from point to point - but there is no chronology given!
The "ID"s seams to be chaotic and there is no date/time information.
How do I get lines from that layer?


texts: ERROR001409 failed to write the output location

texts:
When I try to execute "Summarize Within" tool, it gives me the "001409" error ("Failed to write the output location").

What is the reason of this error and how I can avoid it?


texts: How to execute query based on previous query result and a field filter in arcgis flex api?

texts:
I'm trying to execute recursive queries and each one depend on the previous result. I'm looking for a method to use a filter field in the query execution based on the previous query result.
Is this possible in arcgis flex api? If yes then what method is used in Flex?

here is My DoQuery function:

            private function doQuery(geom:Geometry):void
            {
                try
                {
                    var query:Query = new Query();
                    query.spatialRelationship = Query.SPATIAL_REL_INTERSECTS;
                    query.geometry = geom;
                    query.returnGeometry = true;
                    query.outFields = ['HUC_14'];
                    queryTask.execute(query, new AsyncResponder(onResult, onFault));
                    var GlobalFeaturesArray:Array = new Array();
                    var count:Number = 0;
                    function onResult(featureSet:FeatureSet, token:Object = null):void
                    {

                        if (featureSet.features.length > 0)
                        {

                            count = featureSet.features.length;
                            GlobalFeaturesArray = GlobalFeaturesArray.concat(featureSet.features);
                            var attributes:ArrayCollection = new ArrayCollection(featureSet.attributes);
                            var Collection:Array = new Array();
//In this if clause, I'm trying to append all HUC_14 values from the query to Collection Array in order to run next query.
                            if (count > 0){
                                var i:Number;
                                for (i=0 ; i < count; i++){
                                    Collection.concat(featureSet.features[i].attributes["HUC_14"]);
                                }
                            }

                            while (count > 0 )
                            {
                                var query2:Query = new Query();
                                query2.outFields = ['HUC_14'];
                                Alert.show(Collection.toString(), 'Alert Box', mx.controls.Alert.OK);
                                query2.where = "HU_14_DS in (" + Collection.toString() + ')';
                                queryTask.execute(query2, new AsyncResponder(onResult2, onFault2));
                                function onResult2(featureSet2:FeatureSet, token2:Object = null):void
                                {
                                    count = featureSet2.features.length;
                                    Collection.concat(featureSet2.attributes['HUC_14']);
                                    GlobalFeaturesArray = GlobalFeaturesArray.concat(featureSet2.features);
                                }

                                function onFault2(info2:Object, token2:Object = null):void
                                {
                                    Alert.show(info2.toString(), "Query2 Problem");
                                }

                            }

                            for each (var myGraphic:Graphic in GlobalFeaturesArray)
                            {
                                myGraphicsLayer.add(myGraphic);
                            }
                        }
                        else
                        {
                            Alert.show("No parcels were found", "Try something else");
                        }
                    }
                    function onFault(info:Object, token:Object = null):void
                    {
                        Alert.show(info.faultString + "\n\n" + info.faultDetail, "queryTask fault " + info.faultCode);
                    }
                }
                catch (error:Error)
                {
                    Alert.show(error.toString(), "myGeometryService_simplifyCompleteHandler error");
                }
            }


texts: Combining PostGIS, GeoDjango and OpenLayers

texts:
I'm attempting to create a GIS application to display some data onto a OSM map.

I followed the GeoDjango tutorial to get shapefile WorldBorders data into my Postgis Database, but I am unsure on how I could pull the relevant data into a template for use with OpenLayers.

Any advice on how I may be able to connect this data with OpenLayers?  I've looked around but can't find anything concrete on how to do this, asides from using GeoServer but I would rather stick with my current stack if possible.


texts:
I am using the Gdal-Raster Calculator in QGIS to clip a Input-Raster (3-Bands RGB (A,B,C)) by another Raster (D) (containing 0/1-Values) by using the expression (AD)+(BD)+(C*D), which works well. However the Output contains only one single Band (grey1).
Is there a possibility to get an output with all three input Bands again?


texts: How to silence QGIS 3 nightly build's debug output with PyQGIS?

texts:
When I run the following file sample.py it prints more than 100 lines of debug output. When porting a PyQGIS script to QGIS 3 it is hard to find the expected log statements in between all debug messages.

import os
from qgis.core import *

os.environ["QT_QPA_PLATFORM"] = "offscreen"

QgsApplication.setPrefixPath("/usr", False)
qgs = QgsApplication([], False)
qgs.initQgis()

print(qgs.qgisSettingsDirPath())

qgs.exitQgis()

(Run with python3 sample.py on a nightly installation.)

How can the output be made more readable and the debugging silenced?


texts: 'TypeError: 'tuple' object does not support item assignment' when iterating through the rows of a layer

texts:
I'm trying to assign a value to the attributes of each feature of a feature layer. To do so I am just iterating through each row using 'arcpy.da.SearchCursor'.

#Import geoprocessing.
import arcpy

#Import workspace.
arcpy.env.workspace = r'D:\ArcGIS_Pro_Projects\Optimity_Alf'

#Input parameters
buffer_distance = arcpy.GetParameterAsText(0)

#Input variables.
fc_Opt_All = r'D:\ArcGIS_Pro_Projects\Optimity_Alf\Optimity_Equipment.gdb\Optimity_All'
fc_Opt_Opp = r'D:\ArcGIS_Pro_Projects\Optimity_Alf\Optimity_Equipment.gdb\Optimity_Prospect'
fc_London_Bldngs = r'D:\ArcGIS_Pro_Projects\Optimity_Alf\London_Building.gdb\London_Buildings_Multipatch'
DSM_model = r'D:\Data4Optimimty\ESRI_UK_Bluesky_DSM\ESRI_UK_Bluesky_DSM.gdb\DSM_25cm_2015'

arcpy.MakeFeatureLayer_management(fc_Opt_Opp, 'fl_Opt_Opp')
field_names = ['LOCATION', 'ESTATE', 'ADDRESS', 'OFFSETA', 'OFFSETB', 'Z', 'OFFSETC', 'AZIMUTH1', 'AZIMUTH2', 'VERT1', 'VERT2', 'RADIUS1', 'RADIUS2', 'STATUS', 'CONNECTION']
observers = []
targets = []
#Calculating LOS for each point
with arcpy.da.SearchCursor(fc_Opt_Opp,field_names) as cursor:
    for row in cursor:

        #Adding default values to observers' attributes
        row[1] = 'Prospect'
        row[3] = 1
        row[4] = 0
        arcpy.AddSurfaceInformation_3d('fl_Opt_Opp', DSM_model, 'Z', 'BILINEAR')
        row[6] = row[4] + row[5]
        row[7] = 0
        row[8] = 360
        row[9] = 1
        row[10] = -16
        row[11] = 0
        row[12] = 2000
        cursor.updateRow(row)

However, I am getting the following error:

Traceback (most recent call last):
  File "D:\ArcGIS_Pro_Projects\Optimity_Alf\Arcpy_Scripts\LOS_Many_Points.py", line 28, in <module>
    row[1] = 'Prospect'
TypeError: 'tuple' object does not support item assignment
Failed to execute (LOSManyPoints).

I really can't see what tuple is the error referring to. Besides, if I do the same thing but for just one feature (not iterating) I do not get the error and it works fine.