texts: Using PyQGIS to set a mem_layer values to 2 decimal points

texts:
I am creating 12 mem_layers through a for loop and then using the zonal statistics tool for each one, to calculate a sum value. This sum value has many decimal places and I would like to reduce the amount of decimal places for each value, so that it contains two decimal places. This is the code I am using the create and apply zonal stats to the mem_layers:
# Creating the distance buffers
layer = iface.activeLayer()
feats = [ feat for feat in layer.getFeatures() ]
epsg = layer.crs().postgisSrid()
uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer&field=x:real&    field=y:real&field=point_id:integer""&index=yes"

colors = ["#FCFF00", "#FFF700", "#FFE000", "#FFCC00", "#FFB900", "#FF9E00", "#FF9200", "#FF7B00", "#FF6C00", "#FF4500", "#FF0000", "#FF0000"]

#######
# this list will save all mem layers
mem_layers = []
for color_idx, (intensity, dist) in enumerate(DistInt.items()):

    if dist <= 0:
        continue

    mem_layer = QgsVectorLayer(uri,
                           f'MMI {intensity} Buffer',
                           'memory')
    prov = mem_layer.dataProvider()

#Using zonal statistics to work out affected population sum column
    prefix = rlayer.name()
    params = {'INPUT_RASTER': rlayer.name(), 'RASTER_BAND': 1, 'INPUT_VECTOR': f'MMI {intensity} Buffer', 'COLUMN_PREFIX': prefix+'_', 'STATS': 1}
    processing.run("qgis:zonalstatistics", params)


texts: Plotting Vector Arrows from U and V Component with ArcGIS for Desktop?

texts:
I have two raster layers representing the U and V components of a vector wind. I could generate the magnitude raster by applying Square root in Map Algebra in ArcGIS.

I also generated a raster layer of wind direction by two formulas:

x = ATan2(-v, -u) / (2*3.14159265) + 1

(x - Int(x)) * 360

How can I plot Vector Arrows from U and V Component with ArcGIS for Desktop?


texts: QGIS 3 calculate totals and display on the map

texts:
I have recently started using QGIS 3, at a very basic level really, I am displaying radio coverage drive test logs which are in tab delimited format.
I have sucessfully imported the file, and written several expressions saved into a 'style' to display the data on a base map with a different coloured spot for each signal level.

Within the log file are entries for 'call start', 'call end', and 'call drop, I have also written expressions to display each of these on the map as a different symbol.
If I click on the log file I imported, and select 'show feature count' I am able to see how many calls were started, how many were completed and how many dropped.

What I now want to do is take these totals, perform simple maths on them to create a percentage of successful calls and then display that onto the map alongside the drive test coverage log, in like a title block if you see what I mean.

Would anyone be able to give me some idea of how I would go about that please?
I have tried adding a new scratch layer and then I was hoping to write a similar expression to what I already have but to count the number of call start and call end entries but I cannot see how I reference the data from the previously imported coverage log layer .

This is what the imported log file looks like.

Timestamp   Longitude   Latitude    Level   Qual    EVENT
2018.04.06_09.04.13 -1.138316   52.055508   -92 -10 LOG START
2018.04.06_09.04.13 -1.138316   52.055508   -92 -10 CALL_END
2018.04.06_09.04.13 -1.138316   52.055508   -92 -9
2018.04.06_09.04.15 -1.138316   52.055508   -92 -9  CALL_START
2018.04.06_09.04.22 -1.138316   52.055508   -92 -9  CALL_END
2018.04.06_09.04.24 -1.138316   52.055508   -91 -8  CALL_START
2018.04.06_09.04.31 -1.138316   52.055508   -87 -10 CALL_END
2018.04.06_09.04.34 -1.138316   52.055508   -85 -9  CALL_START
2018.04.06_09.04.41 -1.138316   52.055508   -85 -9  CALL_END

Apologies for the formatting, but each line represents an 'event' and the GPS location and the event which can be call start, call end, call dropped or call blocked.
In the example above, I was not moving which is why the Lan / Lon are all the same.


texts: I am trying to convert from dxf to shp. Any text layer when converted has attributes including text, height, and angle. The angle is in radians, and I can see that it is in degrees in the dxf file. QGIS displays the text as if the angle was in degrees. How can I fix this?


texts: Why does QGIS shift the raster during raster calculator operation

texts:
I have a raster with bathymetry data and need to add 0.6 value to each cell. The operation itself is simple - I just use raster calculator to add the value using following expression: bathymetry@1+0.6. The calculation runs successfully and at a first glance, everything is OK.

However, when I looked closely, it seems that the raster has been shifted on x and y axis, by approximately 0.2m, this value however varies.It looks like this (green is the raster raised by 0.6, gray is the source one):

So the question is, why does this happen? And also, how to avoid it?

@EDIT: I've noticed that this happens for any rasters that are in WGS84 CRS (EPSG 4326). When I carry out the same operation on layers in local, meter-based coordinate systems, everything is fine.

I'm using QGIS 3.14.3 and Windows 10.


texts: Creating custom Coordinate System in PostGIS

texts:
I have a Wayne county Michigan parcel dataset with a custom SRID:
when I bring it into ArcGIS
Projected Coordinate System:    NAD_1983_HARN_StatePlane_Michigan_South_FIPS_2113_Feet_Intl
Projection: Lambert_Conformal_Conic
False_Easting:  13123359.58005249
False_Northing: 0.00000000
Central_Meridian:   -84.36666667
Standard_Parallel_1:    42.10000000
Standard_Parallel_2:    43.66666667
Latitude_Of_Origin: 41.50000000
Linear Unit:    Foot

Geographic Coordinate System:   GCS_North_American_1983_HARN
Datum:  D_North_American_1983_HARN
Prime Meridian:     Greenwich
Angular Unit:   Degree

I am trying to insert this SRID into PostGIS. Is this enough info to go and create the insert statement into PostgreSQL?


texts: PyQGIS - Save processing output to temporary file which is discarded when QGIS is closed (on Windows)

texts:
QGIS 3.10.1, Windows 10:

When using the GDAL tool Clip Vector by Extent from the processing toolbox and saving the output to a temporary file, the output layer appears in the layer list with a special icon (looks like RAM) and the information: "Temporary layer only! Contents will be discarded after closing QGIS":

Is it possible to create such a tempory layer also on Windows when using the tool in PyQGIS?

On Linux this is possible, just save the output to /tmp/.

However on Windows, saving the output to the Windows Temp folder (.../AppData/local/Temp/), the created output is not discarded after closing QGIS.

Here is the code I'm using in PyQGIS:

clipped_layer = processing.run("gdal:clipvectorbyextent", {'INPUT': id, 'EXTENT': canvasExtent, 'OUTPUT': outputFolder + 'Canvas_Clipper.A_' + map.name() + '.gpkg'})

newLayer = QgsVectorLayer(clipped_layer['OUTPUT'], 'A_' + map.name(), 'ogr')
QgsProject.instance().addMapLayer(newLayer, False)
clippedGroup.insertLayer(1, newLayer)

Other answers here on gis.stackexchange are suggesting to use NONE or 'memory:' as output parameter, but this gives an error in PyQGIS. Seems like for GDAL tools there is a need to specify an output path.

When using the tool from the processing toolbox in QGIS, the output is saved to a folder which is prefixed processing_ (see the following screenshot). This folder is discarded when QGIS is closed and I assume this folder is created when starting QGIS (the name of the folder changes for every QGIS session).

Is there a way to use/get the name of this tempory processing_ folder as output folder from PyQGIS when running the GDAL tool or to create such a folder which is automatically discarded when QGIS is closed?


texts: How to achieve proper alpha blending using GDAL command line tools only

texts:
With only using GDAL command line tools is it possible to avoid those white borders around edges of features?

I am using gdal_translate to convert vector PDF to raster mbtiles (OGR is not recognizing PDF format but GDAL does)

gdal_translate -of mbtiles --config GDAL_PDF_DPI 300 test.pdf test.mbtiles
gdalbuildvrt -srcnodata 255 output.vrt test.mbtiles

(I am using first gdal_translate with some -gcp points to georeference nonspatial-pdf but I have removed it here for easy reading)

I am getting this result:

but this is what I want:


texts: How can I translate a GML definition of SRS to a usual format (e.g. WKT or PROJ.4)?

texts:
From GDAL 1.9.0 a new utility is available: gdalsrsinfo.

Although it can yield an XML (GML-based) definition from a source equivalent format or file, it seems that it can only read usual GDAL/OGR forms (complete WKT, PROJ.4, EPSG:n or a file containing the SRS), so no GML definitions.
Indeed:

$ gdalsrsinfo -o xml '+proj=longlat +a=3396190 +b=3396190 +no_defs'
<gml:GeographicCRS...

but:

$ gdalsrsinfo -o wkt '<gml:GeographicCRS...</gml:GeographicCRS>'
ERROR 1: ERROR - failed to load SRS definition from <gml:GeographicCRS...

Any of you is aware of a way I could make this translation?

This question is somehow linked to: Can GeoTools parse GML CRS definitions?


texts: How to create polygon "stamps"  and extend lines by fixed distance? in QGIS

texts:
This is related to Does QGIS have these CAD editing functions?

Cadtools have solved a lot of the issues I had but I need the ability to see how long a line is while I am drawing an orthoganal polygon? Is there a way to do this? can we change the Coordinate output on the bottom toolbar to show meters from last vertex rather than the coordinate or like in ESRI set an extend by X meters at Y angle?

Also can I pre-create a block of say 3m x 1.5m and then drop paste it in into a shape file? Obviously I can copy paste but can you create a stamp so that it fits into angular lines as well...

The current usage is I am trying to create a qgis based system where non-GIS people can create market stalls within a building easily.


texts: Adding vertices to multiple overlapping line segments in QGIS?

texts:
I have a polyline dataset of transit routes for my city. Around the transit stations, up to hundreds of routes overlap completely, but do not follow the correct path that the busses actually travel (possibly due to gps errors). I am trying to correct this and topologically align these segments to a city road network. the dataset in it's raw form is too inaccurate to automatically align via a tool such as the Integrate tool in ArcGIS Desktop.

Attempts at moving the vertices of multiple features simultaneously in ArcGIS Desktop using map topologies or the align features to shape tool resulted in the creation of multipart features, altered the direction of the lines, and in general, ruined the geometry of the lines. the digitization direction of each route is important for routing applications with the dataset.

In QGIS I could simply select overlapping vertices and move them together. This did not produce multipart features and preserved direction. However, the lines in their raw form do not have enough vertices to accurately be aligned to the road network.

Is there a way I can approach this, such as adding vertices to multiple overlapping polylines at the same time in QGIS?

I cannot fully automate this job, as the raw data is too misaligned from the road network.


texts: Finding users near the google search location result

texts:
I have a leaflet map. In the map i have a google search where you type a location name and then the map zooms and put a marker according to the location.
I also have a database table with user ids and their lat/lng .

Now what i want to do is to set in a variable the lat/lng of the marker of google search. Then do a query at my database and get all the users lats/lngs . And if a user is close to the google's marker's lat/lng display him on the map as a custom marker.
(if he is close will be calculated with the distanceTo() function)
i have knowledge of the php query but i stuck in javascript

how do i continue from the google search marker? How i set a variable with his lat/lng?

javascript :

// google location search
 var GoogleSearch = L.Control.extend({
      onAdd: function() {
        var element = document.createElement("input");

        element.id = "searchBox";

        return element;
      }
    });
(new GoogleSearch).addTo(mymap);

    var input = document.getElementById("searchBox");

    var searchBox = new google.maps.places.SearchBox(input);

    searchBox.addListener('places_changed', function() {
      var places = searchBox.getPlaces();

      if (places.length == 0) {
        return;
      }

      var group = L.featureGroup();

      places.forEach(function(place) {

        // Create a marker for each place.
        console.log(places);
        console.log(place.geometry.location.lat() + " / " + place.geometry.location.lng());
        var marker = L.marker([
          place.geometry.location.lat(),
          place.geometry.location.lng()
        ]);
        group.addLayer(marker);
      });

      group.addTo(mymap);
      mymap.fitBounds(group.getBounds());

    });


texts: I want to add photos in my attribute table. I have created a new field but I am not able to select "raster" in the field type. They only show "short integer", "long integer", "float", "double", "text" and "date".


texts: Mapbox Studio - How to slice a rail line shapefile tileset into different parts?

texts:
Mapbox newbie here. I am trying to build a railway map and when its security systems are expected to be completed. So I have added a rail line shapefile as a tileset to a Mapbpx Studio style. I want to create a map where different parts of the line have different colors and when clicked, present a different date. But I am at a loss. I just get a continuous line and that's it.

Wondering how should I go about it since I can not understand how to slice the tileset, which is a line, into different sections, and add data to it.


texts:
I've done tons of work within ArcGIS to analyze and process rasters. I'm pretty well knowledgeable on that subject. A question has come in from my coworker who is attempting to populate elevation data in his SQL server. The first question he asked was if its possible through code to connect to a raster and query it for data, kind of like connecting to an ODBC data source? The idea would be to query for data in an extent and then create records of elevation points within that area. Is that possible with a data driver?

Second, if he has a current version of SQL Server I believe that has some capability to store GIS data correct? Is there a possibility to load the raster into the SQL server as a table and then query it as you would any other table? Or does he require some GIS framework to do this?

Abel


texts:
I want to reproduce on a GIS like Mango, QGIS, Kibana ... the geodata to reproduce the map of this website. I don't know if it is a simple GeoJSON or

I tried a modification to oktogen's answer :

from bs4 import BeautifulSoup
import requests
import json

url = 'https://data.electionsportal.ge/en/event_type/1/event/38/shape/69898/shape_type/1?data_type=official'
r = requests.get(url)
soup = BeautifulSoup(r.text)
splits = [item.split('=',1)[-1] for item in str(soup.script).split(';')]
filtered_splits = [item.replace('"','') for item in splits if 'json' in item and not 'xxx' in item]
links_to_jsons = ["https://data.electionsportal.ge" + item for item in    filtered_splits]
for item in links_to_jsons:
    r = requests.get(item)
    geojson = r.json()
    print(url.rsplit('/', 1)[-1])
    with open(item.rsplit('/', 1)[-1] +'.geojson', 'w') as f:
        json.dump(geojson, f)

And it looks like it was at least one of the two files and it looks like GeoJSON :

{'type': 'FeatureCollection', 'features': [{'type': 'Feature', 'geometry': {'type': 'Polygon', 'coordinates':
[[[4854546.154083, 5284008.666675], [4855951.745246, 5283252.476637], [4856871.642462, 5283565.511518], [4857608.037316, 5283222.205473], ...
..."properties": {"id": 69898, "parent_id": null, "common_id": "Georgia", "common_name": "Georgia", "has_children": true, "shape_type_id": 1, "shape_type_name": "Country", "value": "No Data", "color": null, "title_location": "Country: Georgia"}}]}

But when importing it to Kibana I have the following layer:

I think the problem is a projection error. Kibana must wait to give the project on a map of the world and these were made for a map of borders in the Caucasus.

On Mapbox I have the following error:

Upload errored a few seconds ago
Close
1-8q9493
bounds west value must be between -360 and 360

I don't know how to reset the projection on a GeoJSON file.

However, in QGis everything looks perfect:

So I know I can reproject from QGis ...


texts: Drawing rows in field at the same distance from each other

texts:
I have to recreate on QGIS what you can see in the picture below. I have to draw those lines (tree lines) in my field. The rows need to be at exactly 10m far from each other and oriented North-South. I've already created the polygon of the field, but I don't know how to go further.
Could you help me with this?


texts:
I am writing a plugin in QGIS with python where the user can select one or more objects to view information on one of these objects. I am making it so that the user can click them back and forth between the selected objects. The problem is that you can't see in the map which object you are currently looking at.

Is there a way to mark one of the selected object more or in another way than they already are?


texts: QGIS RGB image to HSL (Hue Saturation Luminance)

texts:
is it possible using QGIS (2.10) on windows to convert an RGB image file to HSL (Hue Saturation Luminance) format.

I've seen various formulas for this written in Java however they only seem to be for converting single values. I am looking to convert whole images.

I also have access to ArcView 10 + spatial analysis and Global Mapper 16.2.


texts: How to fix spatial reference error on CopyFeatures ArcGIS tool?

texts:
When I run the CopyFeatures tool to copy features from an ArcSDE geodatabase to another empty ArcSDE geodatabase (no features, no tables), I get the following result:

The error message is about the spatial reference of the output feature class. The output is created by the tool itself. Why wouldn't it set the correct spatial reference?

Do you have any hint on how fixing this?


texts: Seeking Javascript API for weather data?

texts:
I'm looking for a good javascript api which gives access to a weather data-base.
I want to extract information like average number of rain days, amount of sunshine and average temperature.

Is there an api or an openly available database with information like that available?


texts: How to read Binary overlay (*.ovl files) data programmatically

texts:
Currenlty, Im dealing with binary Overlay files (*.ovl files) and I needed to read them byte by byte and display the content in a possible way.

I came to know that there can be ASCII and Binary format of Overlay content (*.ovl files) and there are few tools which can read ASCII Ovl files. But, could not get anything on Binary Overlay files.

Could you please anyone point me to possible way of displaying the binary overlay content either by any tool or by programatically (Final aim to program it).


texts:
I tried to install QConsolidate by placing the plugin folder ("qconsolidate-master") in:
c:\users\.qgis2\python\plugins\

but received an error "Couldn't load plugin qconsolidate-master due an error when calling its classFactory() method" when starting QGIS. I am using QGIS 2.4.0-Chugiak.

Is there an incompatibility issue?


texts:
I used Pathfinder Office 5.60 to export SSF to RINEX. But the RINEX contains only L1 code and carrier data. It does NOT seem to export L2 code and carrier. BTW I use a GeoXH 6000 that also records GLONASS.

Is there a way to get the whole lot of raw L1 and L2 data from SSF into RINEX, so I can post-proccess it in another software?


texts: Where are Object class extension GUIDs stored in Oracle?

texts: Where are the GUIDs for object class and workspace extensions stored in SDE for Oracle? I've looked through GDB.ITEMS but can't seem to find it. I think SQL Server stores it in an .xml column but I don't see the equivalent in Oracle.


texts: Lines not rendering correctly in ArcGIS Online

texts:
I have some contour lines in a geodatabase. They are z enabled. I'm using ArcGIS Pro to create a scene and then sharing that to ArcGIS Online. I'm using the default elevation layer and have used Pro to project all contour line files to Web Mercator prior to sharing. The symbology of the lines is a simple line symbol. They look great in ArcGIS Pro:

They look awful and are effectively unuseable in ArcGIS online:

So far I've tried altering the line width keeping in mind desired resolution as described in the ESRI docs here.

I've had this problem previously when exporting from ArcMap to PDF and the ESRI docs mentioned something about windows graphics and that Pro should be used instead (can't find the link now). This worked for exporting to PDF. The same is not true of exporting to ArcGIS Online.

Pro doesn't seem to have any options to choose resolution or change anything much when sharing to ArcGIS online and I thought perhaps the coordinates are being rounded? Not sure but would very much like to hear any thoughts or opinion. Geonet has yielded no results so far.

Update

I've now tried sharing the map as a feature service from ArcMap. Viewing it in map viewer on ArcGIS Online the lines are fine. Viewing the same feature service in Scene Viewer in ArcGIS online the same problem exists.

I'd also be interested to know if anyone can replicate this problem.


texts: Symbology Method. "Size" not showing up in "Method", QGIS 2.10.1

texts:
QGIS 2.10 is supposed to be able to allow you to use size instead of color for a gradient method, as seen here: http://anitagraser.com/2015/06/21/qgis-2-10-symbology-feature-preview/.
However, I get no such option. The "Method" drop down only contains color. Do I need to activate this option somewhere else? I am in QGIS 2.10.1.


texts: Geoserver ImageMosaic index is not updating with cURL

texts:
I am successfully adding a new ImageMosaic of GeoTiffs with the Geoserver REST API

curl -v -u admin:geoserver -XPUT -H "Content-type: text/plain" -d "file://D:/tifdirectory" "http://localhost:8088/geoserver/rest/workspaces/wsp/coveragestores/tifdirectory/external.imagemosaic"

The .SHP index is then created in my "tifdirectory" along with the "tifdirectory.properties" file and everything works well after I setup the layer parameters directly in Geoserver.

My GeoTIFFs are updated daily, so I need to update the footprint index. I tried with this command:

curl -v -u admin:geoserver -XPOST -H "Content-type: text/plain" -d "file://D:/tifdirectory" "http://localhost:8088/geoserver/rest/workspaces/wsp/coveragestores/tifdirectory/external.imagemosaic"

At first, it seems to work since there is no error in the response, I receive:

* upload completely sent off: 72 out of 72 bytes
< HTTP/1.1 202 Accepted

But my ImageMosaic index (.SHP) is not updated. What am I missing here? Is there another way to update to mosaic index?

I think that I could make it work by deleting the layer and coveragestore daily before creating it back. But, when I create the layer I need to manually edit the parameters in Geoserver UI. How can I create it with cURL with all the parameters (crs, BackgroundValues, SUGGESTED_TILE_SIZE, etc) already set?


texts: Choosing spatial reference for GPS data?

texts:
We received GPS data from the vendor.
I want using ArcGIS software to create XY event table and create a point feature class.

What the spatial reference would be for the GPS data?
Is it something I should ask the vendor about?
Is it from the GPS device?


texts: Authentication Error when connecting to WMTS in DSM

texts:
I am relatively new to QGIS and use version 3.4.11 for a university project of Polish archaeological features.

I am trying to get the "Digital Elevation Model ISOK – Hillshade" WMTS Layer from "https://www.geoportal.gov.pl/uslugi/usluga-przegladania-wmts" working, but when I try to connect via the Data source Manager, it shows the error: "...failed: Host requires authentication". I use "http://mapy.geoportal.gov.pl/wss/service/WMTS/guest/wmts/ISOK_CIEN?SERVICE=WMTS&REQUEST=GetCapabilities" to connect to the server.

It worked very well before, but when I tried to connect today, the error appeared. What could be the problem?


texts: ArcGIS Append tool imports attributes incorrectly

texts:
I'm working on a project where I routinely truncate and load data from a source annotation feature class to a target annotation feature class. My method is simple, and runs as follows:

Delete Rows tool on target FC
Append tool to append all data from new source FC to target FC.

Below is a screenshot of some key attributes from the source feature class:

Now, here is a screenshot of the target feature class, after all of the above data has been appended into it:

You will see that:

AnnotationClassID has inverted attributes. "Insert" is appended where "Main" should be, and vice versa.
FontName is now "Arial", when it should be "Verdana"
FontSize is populated with seemingly random number values.
TextString attribute is appended correctly

This is the only set of feature classes I've come across where a truncate and load executes incorrectly, and I can't figure out why, nor do I know where to begin to troubleshoot it.

I am using ArcGIS 10.3.


texts: QGIS Python Console SQL Server Layer Incorrect Attributes

texts:
I am using the QGIS python console to load a SQL Server layer in QGIS using the follow code:

uri = "MSSQL:{MYSERVER};database={MYDATABASE};tables= MYTABLE};trusted_connection=yes"
lpBatch = QgsVectorLayer(uri_lp, "PREM_LP_Batch", "ogr")
QgsProject.instance().addMapLayer(lpBatch)

It loads fine, and I can browse the layer fine. However, upon attempting to save any changes made to the geometries, I get the following error:

  Layer PREM_INT_Batch: OGR error setting feature 170: Error updating
  feature with FID:170, [22001][Microsoft][SQL Server Native Client
  11.0]String data, right truncation(0)

I have found the data types/lengths of the attributes are incorrect. I believe this is leading to the truncation error, as it thinks the values of the fields are too long.

If I load the layer using the UI, the data types/lengths of the attributes are as they should be, and I can save any changes I make to the geometries fine.

Any idea how I can load the attributes correctly through the python console?


texts: How to store point data (height) for each point in a linestring in openlayers

texts:
I am reading GPX files into OpenLayers. I would like to be able to save the edited route back to GPX including the height and time information that the original file contained (or use SRTM data to fill in new elevations). Is there an easy way to do this or is it going to be create a geometry.point to store each vertex of my linestring?

Thanks in advance,
Jon


texts: CQL filters don't work with WFS layer in OpenLayers

texts:
I used a WFS layer from GeoServer in OpenLayers it's shown on the map, but when I want to use a filter it doesn't work; I don't know where the error is.

This is the code:

<script type="text/javascript">
OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url=";
function init() {
//set up a save strategy
var saveStrategy = new OpenLayers.Strategy.Save();

//set up projections
// World Geodetic System 1984 projection
var WGS84 = new OpenLayers.Projection("EPSG:4326");
// WGS84 Google Mercator projection
var WGS84_google_mercator = new OpenLayers.Projection("EPSG:900913");

//Initialize the map
//creates a new openlayers map in the <div> html element with id="map"
var map = new OpenLayers.Map ("map", {
controls:[
//allows the user pan ability
new OpenLayers.Control.Navigation(),
//displays the pan/zoom tools
new OpenLayers.Control.PanZoom(),
//displays a layer switcher
new OpenLayers.Control.LayerSwitcher()
],
//projection: WGS84_google_mercator,
displayProjection: WGS84
} );
// map extent
var mapextent =
new OpenLayers.Bounds(-7.70912,33.65164,-7.39876,33.52465).transform(WGS84, WGS84_google_mercator);

//base layers
var openstreetmap = new OpenLayers.Layer.OSM();
var google_maps = new OpenLayers.Layer.Google(
"Google Maps", {
numZoomLevels: 20
}
);
var google_satellite = new OpenLayers.Layer.Google(
"Google Satellite", {
type: google.maps.MapTypeId.SATELLITE,
numZoomLevels: 20
}
);

//wfs-t editable overlay
var wfs_layer = new OpenLayers.Layer.Vector("Editable Features", {
strategies: [new OpenLayers.Strategy.BBOX(), saveStrategy],
protocol: new OpenLayers.Protocol.WFS({
version: "1.1.0",
// loading data through localhost url path
url: "http:// localhost/geoserver/wfs",
featureNS :  "http://www.opengeospatial.net/cite",
maxExtent: mapextent,
// layer name
featureType: "patrimoine",
// geometry column name
geometryName: "geometry",
schema: "http:// localhost/geoserver/wfs/DescribeFeatureType?version=1.1.0&;typename=cite:patrimoine"
})
});
map.addLayers([openstreetmap, google_maps, google_satellite, wfs_layer]);
map.zoomToExtent(mapextent);
}

function filter(){
            var filterType ="cql";
            //var filter = document.getElementById('zone').value;
            // by default, reset all filters
            var filterParams = {cql_filter: null};
            if (OpenLayers.String.trim(filter) != "") {
                if (filterType == "cql")
                    filterParams["cql_filter"] = "batiment='Bastion de la sqala'";

// In this line my parametre it' batiment I made a test on the batiment wich has a value "Bastion de la sqala"

            }
            // merge the new filter definitions
            mergeNewParams(filterParams);

        }

        function mergeNewParams(params){
            Region.mergeNewParams(params);

        }

and when I type in console map.layers[2].features it show me an error:

Cannot read property '2' of undefined


texts: Downloading Sentinel-1 VH time-series for polygon

texts:
from sentinelhub import FisRequest, Geometry, CRS
    time_interval = (start_date.strftime("%Y-%m-%d"), end_date.strftime("%Y-%m-%d"))

    fis_request = FisRequest(layer=XXX,
                             geometry_list=[Geometry(feature['geometry'], CRS.WGS84)],
                             time=time_interval,
                             resolution=YYY)

    fis_data = fis_request.get_data()

I would like to average VH data for a polygon using the sentinelhub API. I am using the FisRequest API to do that but I am not sure what to specify for layer and resolution and if that is indeed the correct way to do it.


texts: Windows software with no need for graphic viewer to explore WFS service?

texts:
The city of Haarlem (Netherlands) publishes all its open data through a WFS service. To examine if I can build an interesting application with it I would like to explore the data this service produces.

What I am looking for is software similar to generic SOAP clients, that lets me:

point to the service endpoint
select an operation
input request parameters
perform a request
show the response data

I have tried QGIS, but a fresh install of the current version produces an error about a missing qgis_core.dll. And I am not sure QGIS can do what I want, since I am not interested in visualizing the data at this moment.

I don't need a graphic viewer, but something that shows the actual requests and responses because I want to find out how to communicate with such a service in my own code.


texts:
I'm trying to correct some topology errors in my geodatabase, now focusing on two feature classes: grass and buildings. They must not overlap. However, there are plenty of spaces where the grass was generalized, thus, covering entire buildings. There are no automatic fixes for this errors - unfortunately the subtract function deletes the overlapping part of both features, when I'd only like to delete the grass feature class and preserve the building feature class.

In this way, I'm using the "Clip..." function on the Editor Menu to use the buildings to cut out this portion of the grass. However, there's no way to select more than one building at a time to apply this function.

I cannot use the Erase Function from ArcToolBox, because the same inconsistence can happen also in the other way: generalized buildings covering grass areas. Also, I wouldn't like to have an output feature class, once I just need to change the shape of my features, nothing else.

Is there any way I could use the "Clip..." function from the Editor Menu with multiple polygons, by selecting them all at a time?


texts: Trying to make a histogram with ui.Chart.image.histogram() in GEE

texts:
I'm trying to make a histogram with the difference of two NDVI images from different dates, but I'm having trouble with the ui.Chart.image.histogram() function when I try to perform with my large asset (https://code.earthengine.google.com/?asset=users/felolillorobles/Araucania_ODEPA_2019_Buffer) to obtain the histogram, equivalent to 46881041 pixels in the region.

The Error is the following:

  Error generating chart: Image.reduceRegion: Too many pixels in the region. Found 46881041, but maxPixels allows only 10000000.

The error suggests I set MaxPixels to a limit suitable for my computation, or set BestEffort to true. However, when I do this, I must change my original variable (ee.Image) to a dictionary using a reducer. This throws another error when I try to pass this new reduced variable (ee.Dictionary) as an input to ui.Chart.image.histogram(), because the function expects an ee.Image (like my original variable).

Another way that I've been trying this is with an ee.Reducer.histogram(), which allows me to use my large asset, but then, I don't know how to graph with the values contained in the resulting dictionary. Is possible to pass this dictionary object to an ui.Chart.image.histogram()?

Here is my script, the characters in // are other options that i was trying:

// CALCULO DE DNDVI PARA SENTINEL 2A SR

//Imagen Sentinel 2A en la fecha 1
var IMGSentinel2A= ee.ImageCollection("COPERNICUS/S2")
.filterDate ('2016-10-01', '2017-03-31')
.filterBounds (geometry)
.filterMetadata ("CLOUDY_PIXEL_PERCENTAGE",'Less_Than', 5);
var Sentinel2filtro = ee.Image(IMGSentinel2A.max());
var Sentinel2Clip = Sentinel2filtro.clip (geometry);

//BANDAS FECHA 1
var b2_1 = Sentinel2Clip.select('B2');
var b3_1 = Sentinel2Clip.select('B3');
var b4_1 = Sentinel2Clip.select('B4');
var b5_1 = Sentinel2Clip.select('B5');
var b8_1 = Sentinel2Clip.select('B8');
var b11_1 = Sentinel2Clip.select('B11');
var b12_1 = Sentinel2Clip.select('B12');

//Imagen Sentinel 2A en la fecha 2
var IMGSentinel2A2= ee.ImageCollection("COPERNICUS/S2")
.filterDate ('2018-10-01', '2019-03-31')
.filterBounds (geometry)
.filterMetadata ("CLOUDY_PIXEL_PERCENTAGE",'Less_Than', 5);
// .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 10));
 var Sentinel2filtro2 = ee.Image(IMGSentinel2A2.median());
 var Sentinel2Clip2 = Sentinel2filtro2.clip (geometry);

//BANDAS FECHA 2
var b2_2 = Sentinel2Clip2.select('B2');
var b3_2 = Sentinel2Clip2.select('B3');
var b4_2 = Sentinel2Clip2.select('B4');
var b5_2 = Sentinel2Clip2.select('B5');
var b8_2 = Sentinel2Clip2.select('B8');
var b11_2 = Sentinel2Clip2.select('B11');
var b12_2 = Sentinel2Clip2.select('B12');

//Normalized differences
var normdif1 = Sentinel2Clip.normalizedDifference(['B8', 'B4']);
var normdif2 = Sentinel2Clip2.normalizedDifference(['B8', 'B4']);
var dndvi = normdif1.subtract(normdif2);
print (dndvi)

//   var hist = dndvi.reduceRegion({
//   reducer: ee.Reducer.histogram(),
//   geometry: geometry.geometry(),
//   scale: 10,
//   bestEffort: true,
// });

var hist = dndvi.reduceRegion({
 reducer: ee.Reducer.mean(),
geometry: geometry.geometry(),
scale: 10,
bestEffort: true,
 });

 print(hist)

 // Pre-define some customization options.

var options = {
  title: 'HISTOGRAMA dNDVI',
  fontSize: 20,
  hAxis: {title: 'banda'},
  vAxis: {title: 'Frecuencia'},
  series: {
    0: {color: 'blue'},
}};
// Make the histogram, set the options.
var histogram = ui.Chart.image.histogram(hist, geometry, 30)
    .setSeriesNames(['dNDVI'])
    .setOptions(options);
    // Display the histogram.
print(histogram);


texts:
I try to set style settings in QGIS 2.18. However I get QColor error. Here is my code:

layer = iface.activeLayer()

layer.setDrawingStyle('SingleBandPseudoColor')
stats = layer.dataProvider().bandStatistics(1, QgsRasterBandStats.Min | QgsRasterBandStats.Max)
minVal = int(stats.minimumValue + 0.5)
maxVal = int(stats.maximumValue + 0.5)
mean = (minVal + maxVal) / 2
s1 = str((minVal * 2 + maxVal) / 3)
s2 = str((minVal + maxVal * 2) / 3)
item0 = QgsColorRampShader.ColorRampItem(minVal, QColor(10, 100, 10), str(minVal) + ' - ' + s1)
item1 = QgsColorRampShader.ColorRampItem(mean, QColor( 153, 125, 25), s1 + ' - ' + s2)
item2 = QgsColorRampShader.ColorRampItem(maxVal, QColor(255, 255, 255), s2 + ' - ' + str(maxVal))
fcn = QgsColorRampShader(minVal, maxVal)
fcn.setColorRampType(QgsColorRampShader.INTERPOLATED)
fcn.setColorRampItemList([item0, item1, item2])
layer.renderer().shader().setRasterShaderFunction(fcn)
layer.triggerRepaint()

The error is:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "c:/users/mustaf~1/appdata/local/temp/tmpippgsf.py", line 10, in <module>
    item0 = QgsColorRampShader.ColorRampItem(minVal, QColor(10, 100, 10), str(minVal) + ' - ' + s1)
NameError: name 'QColor' is not defined


texts: Some help installing pgrouting_dd on mac os x

texts:
I hope someone may be able to help me but i've been grapling for hours trying to install pgrouting with Driving Time to no success. I have been able to build pg routing no problem.

However if i trying to cmake with the -DWITH_DD=ON flag set when performing the make I get the following and it fails:

[ 10%] Building C object         extra/driving_distance/src/CMakeFiles/routing_dd.dir/alpha.o
[ 20%] Building CXX object extra/driving_distance/src/CMakeFiles/routing_dd.dir/alpha_drivedist.o
[ 30%] Building CXX object extra/driving_distance/src/CMakeFiles/routing_dd.dir/boost_drivedist.o
[ 40%] Building C object  extra/driving_distance/src/CMakeFiles/routing_dd.dir/drivedist.o
Linking CXX shared library ../../../lib/librouting_dd.dylib

Undefined symbols:....

Does anyone have any ideas on how to overcome this?

I should also point out that I am using os x 10.6.x

C


texts: Grouping village points based on distance and population size using ArcGIS Desktop?

texts:
I want to group a number of villages based on their distance to each other and total cumulative population size? any suggestions?

I have a vector file of say 100 village points attached to atable containing data on population size and so forth. I want to group these villages such that each group contains 5-6 villages with the cumulative population of size of 5000. Hence, distance and cumulative population size should be used in my clustering of villages. I have used grouping in arcgis 10.2 with no success. I have arcgis 10.2.


texts:
I'm trying to interpolate a point file that comes from a surface water model output. Each subgrid point has an elevation, and that was used with the cell water levels to generate values for percent of period record inundated. I want to smooth out the blocky shape for a few reasons that aren't important to get into here. I'm using Qgis for a few things that I'd need a more expensive license in Arc.

So here's my grid of points. Cancel that, I'm only allowed to post two images. The points are 130' apart.

I set it to 25' resolution and ran the interpolation.
This is what resulted.

The right third is beautiful and exactly what I want, but the left two thirds are all screwed up. I thought it might be the non-rectangular shape, so I exported a perfectly rectangular set but that didn't help. I tried adjusting everything I could think of to no avail. I tried to manipulate the projections, but I'm not very good with Q and couldn't figure that out.

The one thing I did that seemed to elucidate it was I cut out a section of points to the left of the part that worked, and tried that.

Amazingly that worked. It looks like the TIN basically only will do a certain width. I may have a lot of these to do and really don't want to have to chop them all into pieces and then glue them back together. There's got to be a solution to this. Any ideas?


texts: DSpat::create.lines and spatstat::rlinegrid: Adding too many lines?

texts:
Does anybody know of an issue with either DSpat's create.lines or spatstat's rlinegrid? Whenever I use these functions on my Windows (owins), the output has far too many lines.

For example, if I run

DSpat::create.lines(owint, nlines=80, width = 1, spacing=NULL,angle=0)

The output has approx. 2000 lines.

If I do a bit of math so I can set the number of lines output from rlinegrid, based on window diameter:

  dot <- diameter.owin(owin)

  d8 <- (dot / 80)

  rlgt <- rlinegrid(angle = 0, spacing = d8, as.owin(owint))

The output also has approx. 2000 lines.

Simplifying the owin by converting to sp and using gSimplify, then running through the above functions works for the most part, but sometimes the functions still output too many lines.

The shape I am adding lines to is similar to this:

Am I misunderstanding how rlinegrid and create.lines work?

Note: The code to import the image and run through all the code previous to this issue is very long and arduous, but I can upload it if an example is needed.


texts: Display selected feature table in layout

texts:
l have table and want to display selected feature in layout (print version). When select other country (look image) show me these country data in print version. Mostly do this by copy and paste in layout. But l want to do this automatic.
Does anyone know how to do this?
Probably need script tool.

l have tried in this link (Adding table of certain attributes to map layout?) to get answer, but l could not identify what mean there.


texts:
I'm trying to use gstat to predict spatially autocorrelated random variables. I'm using these sea surface temperature data. I first fit a variogram model and use this in predict. The approach works fine if I use the global dataset, but fails with some regional subsets.

My code is below:

env.dat.basin <- env.dat
# remove duplicate locations
env.dat.basin <- env.dat.basin[!duplicated(env.dat.basin[,-c(3,4)]),]

# convert to a spatial points df
coordinates(env.dat.basin) <- ~x + y
proj4string(env.dat.basin) <- CRS("+proj=longlat +datum=WGS84")

# fit variogram
v <- variogram(SST~1, env.dat.basin)
fit.v <- fit.variogram(v, model = vgm(psill = 30, model = 'Gau', range = 5000, nugget = 0), fit.sills = T, fit.ranges = T)
plot(v, fit.v)

# define a gstat object (spatial model) based on variogram fit
g.dummy <- gstat(formula=SST~1, data = env.dat.basin, beta = 1, dummy=T, model=fit.v, nmax = 30)

# make nsim simulations based on the stat object
yy <- predict(g.dummy, newdata=env.dat.basin, nsim=1)

this works, and if I subset the global dataset

env.dat.basin <- env.dat[which(env.dat$ocean == 7 | env.dat$ocean == 513), ]

it works fine too, but with a different subset

env.dat.basin <- env.dat[c(which(env.dat$ocean == 11), which(env.dat$ocean == 7 & env.dat$y <10)), ]

which looks like this:

the prediction returns many warnings that cause R to crash if nsim >1

Warning messages:
1: In predict.gstat(g.dummy, newdata = env.dat.basin, nsim = 1) :
  singular simulation covariance matrix
2: In predict.gstat(g.dummy, newdata = env.dat.basin, nsim = 1) :
  Covariance matrix singular at location [6.3,-21.4333,0]: skipping...
3: In predict.gstat(g.dummy, newdata = env.dat.basin, nsim = 1) :
  singular simulation covariance matrix
4: In predict.gstat(g.dummy, newdata = env.dat.basin, nsim = 1) :
  Covariance matrix singular at location [-10.6,-34.45,0]: skipping...
5: In predict.gstat(g.dummy, newdata = env.dat.basin, nsim = 1) :
  Covariance matrix singular at location [-14.885,-15.5333,0]: skipping...

I read elsewhere that this could be due to the presence of duplicate observations, but I think I removed those and the code works for the global dataset.


texts:
I am using QGIS 2.14.9 on windows 10 to view our geospatial data and use the QuickMapServices plugin connected to 'OSM Standard' as a base map. Recently, I had need to use WFS layers, which are not compatible with our organisation's proxy server. If I do not use the proxy server then QuickMapServices will not download the tiles.

To overcome this I am using GDAL TMS driver with this xml:

<GDAL_WMS>
  <Service name="TMS">
    <ServerUrl>http://tile.openstreetmap.org/${z}/${x}/${y}.png</ServerUrl>
  </Service>
  <DataWindow>
    <UpperLeftX>-20037508.34</UpperLeftX>
    <UpperLeftY>20037508.34</UpperLeftY>
    <LowerRightX>20037508.34</LowerRightX>
    <LowerRightY>-20037508.34</LowerRightY>
    <TileLevel>19</TileLevel>
    <TileCountX>1</TileCountX>
    <TileCountY>1</TileCountY>
    <YOrigin>top</YOrigin>
  </DataWindow>
  <Projection>EPSG:3857</Projection>
  <BlockSizeX>256</BlockSizeX>
  <BlockSizeY>256</BlockSizeY>
  <BandsCount>3</BandsCount>
  <Cache />
</GDAL_WMS>

However, the displayed tiles are so 'blocky' that they are very hard to read and looks poor.

I noticed that if I right click on 'OSM Standard' layer in the QGIS Layers Panel, select properties and uncheck 'smoothing' then the QuickMapServices tiles look just as blocky.

Is there a way to 'smooth' the tiles in the GDAL XML or by using QGIS?

Update: CRS is 3857. Example images added:


texts:
I have 2 Polygon layers. 3.5 million records of all NJ Parcels and a 1 record state wide polygon shapefile of NJ. My goal is to emulate the ArcGIS erase analysis too where I would essentially erase all the areas in the NJ shapefile where the 3.5 million parcels intersect and what I would be left with is a statewide polygon of all NJ roadways.

Here is the query I have put together, it seems fairly straightforward but my question is threefold. 1. should I be using st_intersects instead of && 2. is my query constructed correctly so it will achieve my results? 3. is there a faster way for me to construct this query? indexes, etc.. because it has been running for about an hour already.

create table streetbyerase as
select st_difference(parcels.geom,nj.geom) as geom
from parcels,nj
where parcels.geom && nj.geom


texts: Getting distance from nearest stream and elevation to nearest stream?

texts:
I need to know if there are simple commands for this 2 points:

the distance from the nearest stream, D [m], defined as the length of the path that hydrologically connects the location under exam to the nearest element of the reference drainage network;
the elevation to the nearest stream, H [m], computed as the difference between the elevation of the cell under exam and the elevation of the final point of the aboveidentified path.

For the number 1, i tried "cost distance" but i'm not sure.


texts: OpenLayers 6: Converting extent coords to lon/lat

texts:
I have a bbox extent in the form
[-2460660.814556394,-2935181.886150768,2460660.814556394,2935181.886150768].

I want to convert it into lon/lat coords, probably using proj.toLonLat(), but am not sure how I can do it on an entire extent. Hoping I don't have to split it up, transform it, and then reconnect it again.


texts: Docking tables on top of another and viewing both the tables simultaneously in QGIS

texts: How to dock attribute table in QGIS such that both tables can be docked horizontally and viewed at the same time?


texts: On Windows, where's the location of the pg_service.conf file?

texts:
Trying to install QGIS for a few people but they have Windows. The documentation says to save it on the home folder on Mac but completely forgets about the Windows.

File is needed in order to open the project.

Where to put this file?


texts: I am trying to clip a raster (DEM) with a polygon layer in ArcGIS10.1, using Data Management>>Raster>>Clip. Basically the clipped raster has many pixels outside the polygon boundary, mainly for thin polygons (about 1 pixel wide) or ones with sharp corners.  I have different versions of the original raster, all with 25m cell size (but slightly different outline extent), andwhen I clip them I get slightly different results. (I would expect one or two pixels, but this seems like too much). I'm a fairly new GIS user, and I suspect the problem is something with my environment settings, but I've tried with 'snap to raster' both on and off, and the problem persists. I have to calculate zonal statistics on the clipped raster, and the area calculation shows that the clipped area is not the same (or even close sometimes) as the original polygon. Has anyone encountered this problem or have suggestions to fix it?


texts:
We have a data driven pages enabled MXD. Refreshing this mxd in arcmap works and gives me 9000+ pages.
Via python however, mxd.dataDrivenPages.pageCount returns 0. Even after mxd.dataDrivenPages.refresh().
In ArcGIS resources, they say:

Four elements are required for the creation of Data Driven Pages. They are:
Detail Data Frame,
Index Layer,
Name Field,
Sort Field

In ArcMap I use the FILENAME option in the Name Field. In python, I find that mxd.dataDrivenPages.pageNameField is None. But I'm unable to set this property. Any ideas on what to look for to get my pages via python?
Code is:
import arcpy, oss

logger = oss.log()

mxd = arcpy.mapping.MapDocument(MxdFile)
ddp = mxd.dataDrivenPages
logger.info('Checking required params')
logger.info('dataFrame is : "%s"' % mxd.dataDrivenPages.dataFrame)
logger.info('currentpageid is : "%s"' % mxd.dataDrivenPages.currentPageID)
logger.info('index layer is "%s"' % dpp.indexlayer)
logger.info('pagenamefield : "%s"' % mxd.dataDrivenPages.pageNameField)

logger.info('Refreshing index layer')
ddp.refresh()
logger.info('ddp number after refresh : "%s"' % mxd.dataDrivenPages.pageCount)

This results in:
INFO - Checking required params
INFO - dataFrame is : "<geoprocessing Data Frame object object at 0x0000000010CC
6970>"
INFO - currentpageid is : "0"
INFO - index layer is "VISION.SWG_WPA_MAP_ATLAS_TB"
INFO - pagenamefield : "None"
INFO - Refreshing index layer
INFO - ddp number after refresh : "0"


texts: Creating webapp similar to Opportunity Atlas with ArcGIS Online and Web AppBuilder

texts:
If you have not looked at it yet, the Opportunity Atlas provides analysis on 'Which neighborhoods in America offer children the best chance to rise out of poverty?'
I have been trying to create something similar with ArcGIS Online Web AppBuilder but I am unable to figure out how to create a gis web app that can take inputs similar to Opportunity Atlas and update the map accordingly.
Concretely, the opportunity atlas takes input from users and draw visualizations on the map accordingly by processing the data. The user gives an input, say a county, and the associated census results, and the map lays it out. This will require processing at the database end. However, the webapp builder can only display layers that have already been created. Essentially, I want to create and display layers based on user inputs.
Can you help?


texts:
I have x,y coordinates as excel file.
How can I add them into ArcMap?

Or - I have x,y coordinates of 2 points.
How can I quickly find them AND add/draw them on the map in ArcMap?


texts:
I need to make 30x30 and 20x20 grid layouts with co-ordinate points in the corners. These points then need to be relayed to a GPS to be able to set the grids out in person and then survey them. The results are then needed to be uploaded and georeferenced to a basemap.

(Is QGIS able to easily do this and how simple is it to provide .CSV points for the GPS since I am not sure that it takes .GPX data)

I have managed to make a grid and points, I know the points have and ID, is it possible to make this ID start at 1?

Also can I make the ID numbers appear next to the points?


texts: Placing message in geoprocessing progress dialog?

texts:
I am after a way on using arcpy.GetMessage to achieve the following message in the geoprocessing dialog result window?;

The above 4 datasets have been reprojected to: WGS Web Mercator (auxiliary sphere) and saved in S:\GisData

Does anybody know how this can be achieved?

OK, here's my attempt, which isn't currently working;

# Get the spatial reference
 spatialRef = arcpy.Describe(template).spatialReference.name

#Loop through shapfiles in folder and reproject
 for fc in fcList:
     fcspatialRef = arcpy.Describe(fc).spatialReference.name
     if fcspatialRef != spatialRef:
         arcpy.Project_management(fc, outFolder + "\\" + fc, template)
         projCount = int(arcpy.GetCount_management(fc).getOutput(0))
     else:
         arcpy.CopyFeatures_management(fc, outFolder + "\\" + fc)

# Print shapefile Project results
     arcpy.AddMessage(fc)
# Get the count from GetCount's Result object
 msg = "{0} datasets have been reprojected to: {1} and saved in {2}".format(projCount, spatialRef, OutFolder)
 arcpy.AddMessage(msg)


texts: I am looking for polygons of the London Metropolitan Boroughs before 1965 as kml or kmz (I am working with google maps) or shapefiles.


texts:
I have written some code which offset a certain latitude and longitude with a distance in centimeters. What I am trying to do is based on a coordinate and a vector create two new points along the line described by the vector.

I believe I have succeeded in doing this but the results I am getting are a bit odd. I am not very experienced with map coordinates and projections so I am posting here to see if my results 'look' right.

The following link shows an image I have of two lines:

One extends mostly along the longitude and the other mostly along the latitude. Both lines should have a length of 1 meter. However the red line is noticeably longer than the green line.

The lines are located around longitude 6.080695759673799 and latitude 52.501029191616425 (The Netherlands). I know with the projection I am using (Web Mercator I believe) landmasses become stretched along the poles. Is this what I am seeing with the red line? Is the red line stretched because we are closer to the poles than we are at the equator?

I use a combination of two functions to offset the latitude and longitude. I found the code in this forum post: https://sharpmap.codeplex.com/discussions/395118
I've then rewritten the code to Javascript like this:

function toMercator(latLong)
{
    var x = latLong.y * 20037508.34 / 180;
    var y = Math.log(Math.tan((90 + latLong.x) * Math.PI / 360)) / (Math.PI / 180);
    y = y * 20037508.34 / 180;

    return {x: x, y: y};
};

function toLatLong(mercator)
{
    var lon = (mercator.x / 20037508.34) * 180;
    var lat = (mercator.y / 20037508.34) * 180;

    lat = 180/Math.PI * (2 * Math.atan(Math.exp(lat * Math.PI / 180)) - Math.PI / 2);

    return {
        x: lat,
        y: lon
    };
};

To offset the point I do this:

var lat = 2;
var long = 5;

var xOffsetCM = 50;
var yOffsetCM = 50;

var mercator = toMercator({x: lat, y: long});
// Add offset
mercator.x += normalizedVector.x * (50 / 100); // divide offset by 100 to get meters
mercator.y += normalizedVector.y * (50 / 100);

var latLong = toLatLong(mercator);

After the last function call latLong should contain the new coordinate.


texts: If using ogr2ogr to load a PostGIS DB, does the PG string need host and port arguments?

texts:
The ogr2ogr command in this question
How to Import ESRI Geodatabase format .gdb into PostGIS
was missing the host and port arguments, or so I thought.

My answer included something like this

ogr2ogr -f "PostgreSQL" PG:"host=localhost port=5432 dbname=NY user=postgres" NY.gdb

I know my command works (even if it was down-voted) but I wonder if there are cases when the host and port arguments are not required.

So it's possible for my server to be running on the default port, I see. So, how or why would i set up the server to run on the default port? Thanks.


texts: Merging hundreds of vector layers and sum attributes in QGIS?

texts:
I am still at a fairly beginner level with QGIS.

I am trying to construct a map of mean daily rainfall in Africa over the past ten years. I have access to shapefiles for daily rainfall (with identical columns in the attribute table) and I would like to sum up each day for an entire year period. Then, I want to take the mean of the last ten years and create one shapefile.

So far, I've been working with only two specific days with the thought being that if I can combine two, I can combine 365.

I've tried Vector>Data Management Tools>Merge shapefiles to one, but this just creates one larger shapefile and doesn't sum up overlapping areas. I've also tried  Vector>Data Management Tools>Join attributes by location, which gets closer to what I want, but only keeps the geometry of the target vector layer and loses that of the other.

Let me know if there may be a quicker way to go about creating a shapefile of the ten year average other than what I've proposed.

Below are a couple pictures:

Project with two layers to merge (I added the continent outline for a bit of context)

Sample of layer attribute table (I only need to sum the values from the "Contour" column--this is daily rainfall in mm)


texts: How can I index a scene with no data file for a measurement?

texts:
When indexing a scene that does not have a data file for a measurement defined in the product's definition, what is the best practice for indexing the scene such that ODC returns no data for it (more specifically, the data variable should be filled with the nodata value for that measurement for that time)?
Can we specify an empty string ('') for path for the measurement in the dataset document to do this, or something similarly simple?


texts: How to overlay unprojected image in openLayers 3

texts:
I need to overlay an image that does not move with the map layers. I've considered using the ol.control.Control, as the doc specifically states they can be "informational only," as in my case. All the examples of the control are functional (rotate, etc.), but I need a simplified example that is just the image.

Update: here is what I'm playing with (only testing with DIV, will replace with image), but the DIV does not show as expected:

    var forecastLabel = function(opt_options) {
      var options = opt_options || {};

      var $button = $('<div>').css({
            'background': '#f00'
            ,'width': 100
            ,'height': 100
        })
        .text("BUTTON!");

      ol.control.Control.call(this, {
        element: $button[0],
        target: options.target
      });

    };
    ol.inherits(forecastLabel, ol.control.Control);

      var map = new ol.Map({
          controls: ol.control.defaults().extend([
            new forecastLabel(),
          ]),
   ...


texts: Seeking program, plugin, or combination that allows creation of Voronoi diagrams on globe

texts:
I have a little side project which has been stymying me for a bit. As the title suggests, I want to divide the world up into sections determined by how close they are to one of a series of points, ie a Voronoi diagram. I have technically done this in QGIS, but only in the sense that I have divided the projected map. (If the way I phrased that doesn't make sense, this screenshot should clear things up.
That function is useful for dividing smaller areas into Voronoi diagrams, but since I want to divide the whole world by places' proximity to those points, this is obviously not going to cut it. (Hooray for edge cases.) I know it's possible to create Voronoi diagrams on a sphere; I've seen plenty of examples, came across surprisingly many papers on the subject, and even found a StackOverflow page where someone provided a Python algorithm for spherical Voronoi diagrams. (Also another StackExchange user who seems to have the exact same question, but no answers.)
But I haven't found any plugins for QGIS which would let me do that, nor any other functions or plugins for other similar programs. (Or GIS programs that seem to use spherical geometry at all—looking for that got me a lot of azimuthal orthographic projections, but not much else.)
In theory, I suppose I could learn how to make a QGIS plugin and implement that algorithm I found earlier...but like I said, this is just one step in a little project I'm doing for fun. I'd rather estimate the geodesics and construct the Delaunay triangles myself...but before I start down that road, I figured I'd ask some people who know what they're talking about if the wheel's already been invented and I just haven't found it.
Are there any programs/plugins that would let me do this? (Or barring that, ones which would figure out geodesics?)
Obviously I'd prefer a QGIS plugin over another option if I had multiple choices, but I'm willing to install a different GIS program if I had to.


texts: Maps in Google Earth and QGIS do not match?

texts:
So I am plotting a path of vehicle in the water and the map in Google Earth and QGIS dont line up (atleast with either the NetCDF or the CSV file).

In QGIS everything is in WSG 84 EPSG 4326

Is there a difference between QGIS and Google Earth's coordinate systems? Do I need anything coordinate system?

Heres what I am talking about: The point in QGIS should be more to the right in the yellow part

-- Google Earth --

-- QGIS --


texts: Accessing organization AGOL from Jupyter notebook

texts:
I would like to access my content in my university's AGOL portal using Jupyter notebooks running on chrome in the Anaconda environment.  I have installed the ESRI ArcGIS package in that environment.  I have no trouble logging onto AGOL anonymously by:
gis = GIS()

but when I try logging into my school's portal by:
gis = GIS("*domain*.maps.arcgis.com","*myProfileName*","*Password*")

I get an error:

unable to generate token
invalid username/password


texts: ArcGIS 10.1 SDE database schema other than SDE or DBO

texts:
Just wondering if it's possible to apply a schema other than SDE or DBO to SDE feature or raster tables stored in an enterprise database on MS SQL Server 2008 R2.

I have got a new customise schema and want to apply it to a few of tables in the SDE database.

Any idea?


texts: Are latitudes -90 and +90 both valid in WGS84?

texts:
Geometrically, I understand that a sphere can be addressed in latitude and longitude on the ranges [-90, +90] and (-inf, +inf) respectively.
However, I am working on the specification for plus codes, a geographic addressing system predicated on WGS84. So I am seeking a stronger normative reference than my intuition.
Does WGS84 have any limit on the range of allowable latitudes, and are the exact numbers -90 and +90 allowed (as distinct from "numbers approaching -90 and +90")?


texts: Set access permissions of ArcGIS Server Map Service with ArcPy

texts:
I have a number of Map Services I create with an ArcPy script.
After publishing, I have to go to ArcGIS Server Manager to set the access permissions to "Everyone" by clicking on the heads symbol for every Map Service.

I would like to have an automatic solution to set the permissions with ArcPy after the Map Service is published.

Any idea how this could be achieved?


texts: QGIS Why is the size of pie charts assistant different from symbol size assistant?

texts:
I want to plot a pie chart on a map (layer: BT stations_all). The pie chart circle should be scaled, so that its comparable with other values that take similar ranges (i.e. same legend). In the pie chart scaling is achieved in the "layers properties > diagram > size" dialog and should be equal to the sum of the attributes (pic 1, in this case prey_wwg_m2 + nonprey_wwg_m2, purple+white). In a single value, this is achieved via the symbol size assistant ([pic 2). But I find these dialogues quite different and can't seem to match the two  so that they will take a similar size.
Why are these scale assistants so different for pie charts and simple circles? Is Area=surface scaling and how do I achieve ta comparable scale?


texts: Python labeling expression using numeric field not working

texts:
I'm trying to label parcels from a shapefile in the Python expression below in ArcGIS 10.2. All parcels should be labeled with the last 8 characters of their [PARCEL_NUM], and parcels 10 acres or greater should also include their [PROP_ADDRE]. [Acres] is a Double field.

Instead, all parcels get labeled with [PARCEL_NUM] and [PROP_ADDRE]. Using Maplex, land parcel placement. Any ideas on what's going wrong?

In response to comments, no, no errors verifying the expression. Tried switching to checking for [ACRES] < 10.0 first, no difference. Also, as noted above, [ACRES] is a Double data type. No Null values in any of the three fields.

def FindLabel ( [PARCEL_NUM], [PROP_ADDRE], [Acres]   ):
  if [Acres] >= 10.0:
    return [PARCEL_NUM][-8:] + '\n' + [PROP_ADDRE]
  else:
    return [PARCEL_NUM][-8:]


texts: QProgressDialog in QGIS 3.0 not updating

texts:
I am using a QProgressDialog to show the progress of a routine that is running.  I found this very helpful question/answer how-to-show-a-moving-progress-bar-in-qgis, but I can not reproduce successful results using QGIS 3.0.  I copy and paste the code into the python editor and

replace:
 from PyQt4.QtGui import QProgressDialog, QProgressBar
with:
from PyQt5.QtWidgets import QProgressDialog, QProgressBar

also print sum with print (sum) and then run the script.  The dialog displays with the window title ("Progress") and a blank dialog.

Can you provide a code sample of how to implement a progress bar in a dialog using QGIS 3.0?

below is the code I tried:

from PyQt4.QtGui import QProgressDialog, QProgressBar

def progdialog(progress):
    dialog = QProgressDialog()
    dialog.setWindowTitle("Progress")
    dialog.setLabelText("text")
    bar = QProgressBar(dialog)
    bar.setTextVisible(True)
    bar.setValue(progress)
    dialog.setBar(bar)
    dialog.setMinimumWidth(300)
    dialog.show()
    return dialog, bar

def calc(x, y):
    dialog, bar = progdialog(0)
    bar.setValue(0)
    bar.setMaximum(100)
    sum = 0
    progress = 0
    for i in range(x):
        for j in range(y):
            k = i + j
            sum += k
        i += 1
        progress = (float(i) / float(x)) * 100
        bar.setValue(progress)
    print sum

calc(10000, 2000)


texts: Improving the code to display also the values 0 or null

texts:
The following code calculates the median in the grouped data. I run this code on my feature class with 18,000 features, but the Python IDLE only prints 10,000 results, that is, skips features for which the median cannot be calculated. I would like to improve this code so that Python IDLE also prints results for features for which the median cannot be calculated. For example, to print 0 or null values for these features.

import arcpy
input_table = r'put the path to your table here'
#lcb is the lower class boundary
width_age_range_1 = 14.5
lcb_age_range_1 = 0.0
width_age_range_2 = 50.0
lcb_age_range_2 = 14.5
width_age_range_3 = 35.5
lcb_age_range_3 = 64.5
median_list = [ ]
with arcpy.da.SearchCursor(input_table, ["TOT", "TOT_0_14",  "TOT_15_64", "TOT_65"]) as cursor :

for row in cursor :

    #if all numeric fields hopefully
    n_over2 = (row[0]/2.0)

    if row[1] >= n_over2 :
        #cumulative freq of groups before median group
        cumu_freq = 0.0
        if row[1] != 0.0 :
            median = lcb_age_range_1 + ( ((n_over2 - cumu_freq)/row[1]) * width_age_range_1)
            median_list.append(median)
            continue
        continue

    if (row[1] + row[2]) >= n_over2 :
        #cumulative freq of groups before median group
        cumu_freq = row[1]
        if row[2] != 0.0 :
            median = lcb_age_range_2 + ( ((n_over2 - cumu_freq)/row[2]) * width_age_range_2)
            median_list.append(median)
            continue
        continue

    if (row[1] + row[2] + row[3]) >= n_over2 :
        #cumulative freq of groups before median group
        cumu_freq = row[1] + row[2]
        if row[3] != 0.0 :
            median = lcb_age_range_3 + ( ((n_over2 - cumu_freq)/row[3]) * width_age_range_3)
            median_list.append(median)
            continue
        continue

    else:
        median_list.append("error")
for median in median_list :
print(median)
#obviously can add a field and append these values in rather than printing them. look at arcpy.da.UpdateCursors


texts: Labeling vertical lines with Shape Length = 0 in ArcMap?

texts:
I am trying to add labels to my linear feature class. Some of my lines are vertical (i.e. I only see a point representing the start or end node of the line, because they have the same XY coordinates but different Z coordinates).

Unfortunately, it seems ArcMap cannot understand where to place the label (e.g. along the line, before the start, etc.) because the line itself is null (i.e. has 0 (zero) shape length).

The same seems to happen in ArcGIS Pro as well (although in the 3D viewer this problem does not occur, see the image below).

Is there a way to solve this?


texts: Using PyQGIS for obtaining data from MS4W WFS link

texts:
I have the following WFS link:
http://118.abc.xyz.ef:90/cgi-bin/mapserv.exe?map=/ms4w/amc/wfs_data.map

I want to read and store the data contained in this link using PyQGIS so that it could be automated. Initially to test it

I opened QGIS v3.14.15-Pi
Went to Layer->Add Layer->Add WFS Layer->New
Put in the name with the URL
In WFS Options->Version I clicked detect which results in 2.0.0 (have tried all the other
options as well)
Set 100 (as well as empty) as the page size
Clicked Ok and Connect
It shows two layers named ms:SNPP_NE and ms:MODIS_NE
I added both layers and as it turns out, only the former shows some data

This confirmed that the WFS link does expose some data
After this I decided to use PyQGIS to run a script to add this WFS layer programmatically. Here is what I tried:
import os
from qgis.core import QgsVectorLayer

path_to_airports_layer = "http://118.abc.xyz.ef:90/cgi-bin/mapserv.exe?map=/ms4w/amc/wfs_data.map"
vlayer = QgsVectorLayer(path_to_airports_layer, "ms:SNPP_NE", "WFS")
#print(vlayer)
#if not vlayer.isValid():
#    print("Layer failed to load!")
#
#else:
    #QgsProject.instance().addVectorLayer(path_to_airports_layer,"ms:SNPP_NE","ogr")

iface.addVectorLayer(path_to_airports_layer,"ms:SNPP_NE","WFS")

#import os
#import urllib
#params = {
#    'service': 'WFS',
#    'version': '2.0.0',
#    'request': 'GetFeature',
#    'typename': 'ms:SNPP_NE',
#    'srsname': "EPSG:3857"
#}
#uri2 = 'http://118.abc.xyz.ef:90/cgi-bin/mapserv.exe?map=/ms4w/amc/wfs_data.map'
#vlayer = QgsVectorLayer(uri2, "my wfs layer", "WFS")
#QgsProject.instance().addMapLayer(vlayer)

Several attempts at trying to add the above WFS layer resulted in the compiler indicating that the layer was invalid (which begs the question of how QGIS was able to load it in the first place). I even tried using the Processing->Graphical Modeler tool in QGIS hoping to generate the PyQGIS script using a model.
In order to do that, under Algorithms, I chose Modeler Tools->Load layer into project, put in the WFS layer URL and clicked ok. But running the model failed as seen below

I even tried to access the GetCapabilities xml/gml document directly from the browser with the following results:

For the above links I referred to the MapServer WFS documentation. What I wish to do is to somehow obtain the data exposed by the WFS service and store it in some form maybe a GeoJSON or even in PostgreSQL. But all my efforts have thus far been fruitless.


texts:
I have a huge CSV file of WGS84 latitudes and longitudes, and I would like to turn these into a large image file that I can print (with each point being a pixel on the image).

What is the best way to do this?

Currently I'm planning to convert them into eastings/northings and then use imagemagick to print a circle at each point, following the approach taken by this author (in the imagemagick script at the bottom of the post).

But I'm wondering if there's a simpler way. I'm not a GIS expert, and perhaps there's even a dedicated geospatial library for doing this?

UPDATE

As requested, here's a sample of the data. The points are all within the UK, but they don't form a regular grid.

Time,latitude,longitude
2014-12-27 18:35:10,51.7438644,-2.2094204
2014-12-27 18:34:08,51.7438657,-2.2094219
2014-12-27 18:33:06,51.7438657,-2.2094219
2014-12-27 18:32:04,51.7438689,-2.2094198
2014-12-27 18:31:03,51.7438716,-2.2094109
2014-12-27 18:30:01,51.7438537,-2.2094211

FURTHER UPDATE

To clarify, I'm looking for a large image that I can print. The points aren't on a regular grid: I was hoping that I could turn each point into a single black pixel on a transparent background.


texts: Where do I find TerraGO GeoPDF hands-on tutorial online?

texts:
Where can I find a TerraGO GeoPDF hands-on online tutorial?

I'd like some examples I can use on my computer to learning how to use it. My agency has their toolbar on our computer system.


texts: Clip shapefile based on another shapefile

texts:
I have two shapefiles with me, one big and one small. I want to clip part of big shapefile where the small shapefile is located. I tried QGIS vector intersect by geometry, but that is returning me empty file.

Is there any command line argument from gdal that can do the job for me?


texts: Presenting Areas of Land from shapefile to 1 decimal position and convert from m-squared to hactares

texts:
I've got a shapefile containing multiple polygons, attached to an attribute table. I'm looking at labeling all these polygons using the attribute 'AREA' column. However, these areas have not been rounded off and are in meters squared whereas I need them to be in Hectares. How would I convert these into hectares then present as a label.

Also, how can I create a white background for these labels?

I'm using ArcMap 10.3.1


texts:
I'm trying to create a map using a basemap pulled in using ggmap (get_stamenmap) with some shapefiles plopped on top (geom_sf). I want to change axis labels! Why's it so difficult!

When I only plot my basemap, I get axes labeled with lat/lon and NO characters - i.e. -19
When I add in the shapefiles using geom_sf I get axes labeled with lat/lon AND characters - i.e. 19°S
When I try change axis labels with scale_x_discrete along with sf_coord(expand = F) I can get what I want
When I add scale_y_discrete with label text as I want them I keep getting this error:

Error: Breaks and labels along y direction are different lengths

What is going on here? Are there invisible breaks on the y axis that I'm missing?

This works:

ggmap(SA) +
  geom_sf(data = traj_outSF, alpha = 0.4, inherit.aes = F) +
  coord_sf(expand = FALSE) +
  xlab(expression(paste("Longitude (", degree,"E)"))) +
  ylab(expression(paste("Latitude (", degree,"S)"))) +
  scale_x_discrete(breaks = c(33.5, 34, 34.5, 35, 35.5), labels = c("33.5", "34", "34.5", "35", "35.5"))

and gets me this:
correct x axis

As soon as I add in scale_y_discrete I get the error message

ggmap(SA) +
  geom_sf(data = traj_outSF, alpha = 0.4, inherit.aes = F) +
  coord_sf(expand = FALSE) +
  xlab(expression(paste("Longitude (", degree,"E)"))) +
  ylab(expression(paste("Latitude (", degree,"S)"))) +
  scale_x_discrete(breaks = c(33.5, 34, 34.5, 35, 35.5),
                   labels = c("33.5", "34", "34.5", "35", "35.5"))
  scale_y_discrete(breaks = c(20, 19.5, 19, 18.5, 18,17.5),
                   labels = c("20","19.5","19","18.5","18","17.5"))

Heres a reprex:

g = st_sfc(st_point(c(34,-19)))
st_crs(g) <- 4326

SA <- get_stamenmap(bbox = c(33.18, -20.3, 35.8, -17.3), maptype = "toner-lite", zoom = 11)

ggmap(SA) +
  geom_sf(data = g, inherit.aes = F) +
  coord_sf(expand = F) +
  scale_x_discrete(breaks = c(33.5, 34, 34.5, 35, 35.5),
                   labels = c("33.5", "34", "34.5", "35", "35.5")) +
  scale_y_discrete(breaks = c(20, 19.5, 19, 18.5, 18,17.5),
                   labels = c("20","19.5","19","18.5","18","17.5"))

Thanks!


texts: Fixed Distance Buffer distance not working on polygons with holes

texts:
I am trying to implement this solution (How to merge single polygons to multi part polygons by maximum distances?), but when I create a buffer around my polygons, those with holes are ignored. I tried to Delete Holes using the processing tool. It does fill the "real" holes but these features:

are not filled in. Polygons with these are also ignored when I try to buffer them...

If I try to fix them manually, I get this:

Also, these polygons were created from a raster using polygonize.
How do I go about being able to create a buffer around such polygons?


texts:
How to query or find the geometry data type a given Oracle Database is using in SQL Developer (i.e. SDO_Geometry, ST_Geometry, etc)?

I thought I could just describe a feature class within the database and get the result, but it just returns the field types. The geometry field just returns 'NUMBER', which isn't helpful.


texts: Error when using gdal/ogr and python to clip a "shp" file by a "shp"

texts:
I am new to the gdal/ogr library for python and I was wondering how to clip a shp file by another shp file. I think I am having trouble with the syntax. Thanks - Stefan

I keep getting a windows error message: "WindowsError: [Error 2] The system cannot find the file specified"

import subprocess

##The features used to clip the input features.**

clipping_shp = "US_States_STATE_NAME__Oregon.shp"

##The feature class to be created.**

output_shp = "output1.shp"

##The features to be clipped.**

input_shp = "broadcast_national_ALL.shp"

##Clipping process**

subprocess.call(["ogr2ogr", "-f", "ESRI Shapefile", "-clipsrc", clipping_shp, output_shp, input_shp])


texts: Leaflet GeoJSON popup marker customization

texts:
I would like to customize leaflet .bindPoPup geoJSON marker.
I found some solutions, although I would like to use my existing one.

Basically the operation leads to merge two functions as per below:

var sitec = {
    "type": "FeatureCollection",
        "features": [
            {
             "type": "Feature",
             "properties": {
                 "Title": "Sitec IS",
                 "Head": "7400 Beach Drive",
                 "Description": "Gavin Sinclair"
        },
            "geometry": {
            "type": "Point",
            "coordinates": [
            0.16964435577392578,
            52.29220753602784
                ]
                  }
                    },
{
  "type": "Feature",
  "properties": {
"Title": "JHG",
    "Head": "Shortstanton sidings",
"Description": "Conor Murphy"
},
  "geometry": {
    "type": "Point",
    "coordinates": [
      0.05458831787109375,
      52.29163006501503
    ]
  }
}
              ]
        };

 //L.geoJSON(sitec).addTo(map);

var geojsonMarkerOptions = {
radius: 8,
fillColor: "#ff7800",
color: "#000",
weight: 1,
opacity: 1,
fillOpacity: 0.8
};

L.geoJSON(sitec, {
pointToLayer: function (feature, latlng) {
    return L.circleMarker(latlng,
 geojsonMarkerOptions);
}
}).addTo(map);

L.geoJSON(sitec, {
    onEachFeature: function (feature,
   latlng, layer) {
        layer.bindPopup('
'+feature.properties.Title+'

Address: '+feature.properties.Head+'

'+feature.properties.Description+'
');
return L.circleMarker(latlng, geojsonMarkerOptions);
}

    }).addTo(map);

If I merge these two functions I got nothing. If I leave them both I have a custom marker, but a .bindPopup remains the same as per in the image below:

Could you help me sort this out?


texts: Is there an alternative to Shapeburst Fill that does not fade out in QGIS?

texts:
I would like to highlight the areas of certain shapes with a thick stroke. However, the accuracy of my shapes is quite high (which is needed for certain zoom details), so that the "simple line" does not look good where more detailed coastal zones or small land areas are present.

The shapeburst fill does a nice job in filling small areas completely, but I don't like the fading effect. Is there a way to use something similar to the shapeburst but without fading?

Alternative thought: Could I add a copy of my shapes on top and colour it with an offset towards the inside in white to create a similar effect?

Image 1: I don't like it.

Image 2: Shapeburst does a good job, but I'd like to have it without fading.


texts:
I will be using a very large dataset (yfcc100m dataset, a dataset containing metadata of Flickr's photos/video) which contains fields longitude and latitude. I'm planning to store the dataset into a MySQL database. Also, I have downloaded shapefiles (.dbf, .prj, .sbn, .sbx, .shp, shx) which contains a country's province borders. I have both shapefiles for the country as a whole and the provinces.

I would like to know if it is possible to query the dataset using the shapefile. What I mean is can I filter out every other data from the dataset which doesn't fall within the borders of the shapefile? If so, can you give me some examples? I'm using Python.

The output is suppose to be a list of all the photos found in the yfcc100m dataset found within the country and grouped by provinces.

The dataset is 14GB while the shapefile is about 9MB. I'm kind of worried since I'm only using my laptop which only has an Intel dual core processor.

I'm new to GIS technology. Please excuse me since I'm not yet familiar into using GIS.


texts:
If I have a PostGIS database that has a table with an integer column that has a sequence that will normally autofill it, how do I get QGIS to not make me fill out that field (and not try to put NULL in it either)?

Everything I seem to do seems to make me either manually manage coming up with new "ID"s or try to put NULL in the column and then fails because the "ID" column is defined as being not NULL.

What do I need to do in QGIS to get it to not specific a value for that column and just let the database (PostGIS) fill it in?

This is with QGIS version 3.8.

This is the definition of one of the tables I'm having trouble with:

CREATE TABLE public.utility_pole (
    id integer NOT NULL,
    pole_owner text,
    pole_owner_primary_label text,
    latlong public.geometry(Point,4326) NOT NULL
);

ALTER TABLE public.utility_pole OWNER TO postgres;

ALTER TABLE public.utility_pole ALTER COLUMN id ADD GENERATED ALWAYS AS IDENTITY (
    SEQUENCE NAME public.utility_pole_id_seq
    START WITH 1
    INCREMENT BY 1
    NO MINVALUE
    NO MAXVALUE
    CACHE 1
);

ALTER TABLE ONLY public.utility_pole
    ADD CONSTRAINT utility_pole_id_pk PRIMARY KEY (id);


texts: Any plans to change GP toolbox format (.tbx file) to text-based format?

texts:
Have you heard if esri has any plans to change GP toolbox format (.tbx file) to text-based format in future?

The issue is that we can't track changes in binary tbx file using source control system. With text based format it would be easy.

Most of our tools are .NET (COM) based. Some used for printing are Python based. Therefore we use a custom toolbox (.tbx), not Python toolbox (.pyt).

Or alternatively (maybe this would be even better approach) esri can provide a tool buildTbx which takes some text-based toolbox descriptor file, finds what tools should be in toolbox, what are inputs and outputs, also other files needed to be included in toolbox (.py files) and  creates that binary tbx file as an output.


texts:
I am trying to return the vertices of a linestring using the sf package. So far I have:

ls  <- st_sfc(st_linestring(rbind(c(0,0),c(0,1))),
        st_linestring(rbind(c(0,0),c(10,0))))

get_line_endings <- function(line){

  res <- lapply(line, function(x){
            pnts <- st_coordinates(x)
            pnts <- pnts[c(1, nrow(pnts)),]
            data.frame(pnts)
          })

  res <- do.call("rbind", res)
  res <- res[!duplicated(res),]

  st_as_sf(res, coords = c("X", "Y"), crs = st_crs(line))
}

pnts <- get_line_endings(ls)

plot(ls)
plot(pnts, add = TRUE)

Is this error prone? Can I rely on st_coordinates being "in order"? Is there a more efficient way of doing this?


texts: Different Leaflet / Folium markers for different zoom levels?

texts:
I'd like to use, say, a circle or map pin when zoomed out, and for this to change to a polygon / shapefile when zoomed in (specifying the zoom level at which the change occurs).
Can Leaflet do this?


texts: Path distance tool only downhill using ArcGIS Desktop

texts:
I'm working on a protection forest model against avalanches and rockfall processes and I have a problem with the path distance tool. I want to calculate the accumulated costs of rockfs travelling DOWNHILL from source cells which represent cliffs.
My cost raster looks like this for not forested terrain and like this for forested terrain:

Slope:          not forested Values:        forested Values:
0 -5 °          10              30
5 - 10°         9               27
10- 15°         8               24
.......         ...             ...
45° and steeper     1           3

If the accumulated path from the start cells is about a certain threshold it will be unlikely that the rockfall process will reach further than this.

How can I use the path distance tool to only calculate the pathdistance for areas downslope from a source cell?

How is it possible to stop the process when it goes uphill again?

For example if process hits valleybottom and opposite slope.


texts: is there a public street database that I can do SQL queries on?

texts: Sorry if this isn't strictly a GIS question. I want to know if there's a public database of US streets and US States that I can simply connect to (e.g. JDBC) and run SQL queries. For example, I want to find all the streets named "California" in the state of California. I've read about Gisgraphy and OpenStreetMaps and Geonames, but I don't see how to do the kind of query I'm talking about.


texts:
I am trying to work out centroid of a polygon, but I need this in British National Grid.
Does anyone know a way of getting it to work in QGIS 3.10?


texts: Rendering CSV file as annotation on QgsMapCanvas

texts:
I am trying to render a CSV file as an annotation layer on QgsmapCanvas using QGIS3 C++ API.
but when I make a vector layer using "delimitedtext" as data provider and test layer->isvalid(). it says invalid.
Here is the piece of code:
QString uri = QString("file:///C:/Users/Puneet/Desktop/annotation/countries.csv?delimiter=,&crs=epsg:4326&xField=latitude&yField=longitude");
qDebug()<< "URI :: " << uri;
QString baseName = QFileInfo(uri).completeBaseName();
qInfo()<< "annotation layer basename::"<< baseName;
QgsVectorLayer *layer = new QgsVectorLayer(uri, "countries", "delimitedtext");
qDebug()<< " annotation layer validity" << layer->isValid();   // returns false.

Here is the image of CSV file.

How can I fix this?  What am I missing?


texts: ST_DWithin Calc Meters - transform or cast?

texts: I want to use ST_DWithin on data which is stored in SRID geometry of 4326 and use meters as the distance parameter. Is it more efficient to do a cast (e.g. data.geom::geography) or a transform to an SRID with units of meters (e.g. ST_Transform(geom, 3857)? Or neither?


texts: ArcGIS 10.1 Legend Issue: Gap when using the "only show classes that are visible in the current map extent"

texts:
I have a problem with a gap that appears when the setting "only show classes that are visible in the current map extent" is in use and when there is a number of items in a row that does not display;

Converting the legend to a graphic and doing it manually is not an option.  If I switch off the Layer files which fall outside of the view extent, the problem goes away, but then it defies the purpose of this handy (yet frustrating) function.

Any suggestions on how to fix this?


texts:
I'm trying to make a GetMap request with a multiple cql_filter as following:

http://localhost:52459/handlerGeoserver.ashx?
SERVICE=WMS&
VERSION=1.3.0&
REQUEST=GetMap&
FORMAT=image%2Fpng&
TRANSPARENT=true&
LAYERS=v_expedientesmapa&
cql_filter=(Year between (2017,2022)) AND IdEstadoPublico IN (1,2,3,4,5,6,7,8,9)&
time=1582632261553&
CRS=EPSG%3A25830&
STYLES=&
WIDTH=1656&
HEIGHT=1361&
BBOX=407077.33718499093%2C4550658.264466829%2C817875.6737958044%2C4888276.945703131

I'm getting filter parsing error. What is wrong with that filter? I'm using the official documentation for this

Edit:

Partially solved adding a double filter (year > 2017) AND (year < 2022)

GeoServer logs:

Caused by: org.geotools.filter.text.cql2.CQLException: Encountered "between ( 2017 ," at line 1, column 7.
Was expecting one of:
    "(" ...
    ":" ...
    <IDENTIFIER> ...
    "." ...
    "*" ...
    "/" ...
    "+" ...
    "-" ...
    ")" ...
    "not" ...
    "in" ...
    "tequals" ...
    "before" ...
    "after" ...
    "during" ...
    "like" ...
    "ilike" ...
    "exists" ...
    "does-not-exist" ...
    "is" ...
    "between" "(" "(" ...
    "between" "(" "[" ...
    "between" "(" <IDENTIFIER> ...
    "between" "(" <DATE> ...
    "between" "(" <DATE_TIME> ...
    "between" "(" "-" ...
    "between" "(" <INTEGER_LITERAL> "*" ...
    "between" "(" <INTEGER_LITERAL> "/" ...
    "between" "(" <INTEGER_LITERAL> "+" ...
    "between" "(" <INTEGER_LITERAL> "-" ...
    "between" "(" <INTEGER_LITERAL> ")" ...
    "between" "[" ...
    "between" <IDENTIFIER> ...
    "between" <DATE> ...
    "between" <DATE_TIME> ...
    "between" "-" ...
    "between" <INTEGER_LITERAL> ...
    "between" <FLOATING_LITERAL> ...
    "between" <STRING_LITERAL> ...
    "between" "true" ...
    "between" "false" ...
    "between" "point" ...
    "between" "linestring" ...
    "between" "polygon" ...
    "between" "multipoint" ...
    "between" "multilinestring" ...
    "between" "multipolygon" ...
    "between" "geometrycollection" ...
    "between" "envelope" ...
    "between" "srid" ...
    "=" ...
    ">" ...
    "<" ...
    ">=" ...
    "<=" ...
    "<>" ...
     Parsing : (Year between (2017,2022)) AND IdEstadoPublico IN (1,2,3,4,5,6,9).
    at org.geotools.filter.text.ecql.ECQLCompiler.compileFilterList(ECQLCompiler.java:122)
    at org.geotools.filter.text.commons.CompilerUtil.parseFilterList(CompilerUtil.java:224)
    at org.geotools.filter.text.ecql.ECQL.toFilterList(ECQL.java:177)
    at org.geotools.filter.text.ecql.ECQL.toFilterList(ECQL.java:159)
    at org.geoserver.util.XCQL.toFilterList(XCQL.java:31)
    at org.geoserver.ows.kvp.CQLFilterKvpParser.parse(CQLFilterKvpParser.java:26)
    ... 102 more


texts: How to avoid slivers when merging (union) two shapefiles in QGIS?

texts:
Here's the kind of tool I want:

Preference A: Union the two polygon shapefiles that where the wetland file takes precedence over the forest file when the distance between the edges is less than 20m.

Preference B: A way to set Eliminate sliver tool to recognize the difference between a narrow sliver and a small round polygon OR to set sliver to only join another polygon with the same attribute in one field.

Situation:

I have two landcover shapefiles (one with forest cover and the other with wetlands) that I want to combine into one single landcover shapefile.  I am running into the problem that the delineations are slightly off so that it ends up with many slivers in the spaces between the delineations.

I start off with the two shapefiles (wetlands in bright green, forests in dull green). It can be seen that they are delineating roughly the same feature here (not always the case, otherwise I could just use one layer) but that they are slightly off. The wetland layer is generally more precise and I would to set it as the primary in cases of conflict/slivering.

This is the result after using the SAGA union tool in QGIS.  (The QGIS union tool available in the Vector menu option creates duplicates).

I've tried sorting the polygons by area covered and then use the Eliminate sliver tool as suggested in this post (How to fill holes in polygons automatically?) but some of my slivers are larger than some of my wetland polygons that I want to keep. Is there a function where the option equivalent to "Merge selection with the neighboring polygon with the" in Eliminate sliver tool is to merge with the neighbouring polygon with the same value in a certain attribute field?


texts:
I am trying to import csv files in ArcGIS to do a join. However, some of the fields are identified as double, some as long and some as long. Therefore the join results in bunch of columns with NULL value. I am not sure why it happens.

I tried to save the numbers with three decimal position in excel as csv format but when I open the number after decimal does not change. Since I have more than 256 columns, I can't save it as xls/xlsx file.

I think using python could be a solution. Any other thoughts? Python experts, could you please provide me any help. Attached is the screenshot..


texts: Set default version of Python in ArcGIS to 64-bit?

texts:
This is the output from Python terminal within ArcMap 10.6.1 for Windows 10:

>>> import sys, numpy, matplotlib
>>> print(sys.version, numpy.__version__, matplotlib.__version__)
('2.7.14 (v2.7.14:84471935ed, Sep 16 2017, 20:19:30) [MSC v.1500 32 bit (Intel)]', '1.9.3', '1.5.2')

How can I set the default Python version to 64-bit? I don't want to change the System Environment Variables as it would probably mess up with the other Python 3.6 from Anaconda3 64-bit installed on the same PC.

Here is my Python27 folder structure


texts:
I am looking for preferably a Linux command line tool that would spit long polylines in my shapefile to shorter segments. Length limit can be defined by the length in meters or in max number of polyline points and doesn't need to be precise.

My problem is that I have a shapefile with contours spanning the whole country and it makes rendering map with Mapnik very slow.


texts: Displaying (or filtering) features in one layer based on features selected in another

texts:
I have points in a layer "Site-Centroids".

I generate a lot of isochrones (polygons) around them, into "Isochrone_Site-Centroids".

With all the isochrones displayed, the map is too busy.

What I want is that no isochrones are initially displayed.

Only when the points from the original "Site-Centroids" are selected, do I want the corresponding isochrones become visible.

The fid's from the "Site-Centroids" are carried over to the isochrones in "ftid", so in the isochrone / polygon layer, I was hoping to use:

is_selected(get_feature('Site-Centroids','fid',ftid))

either as a rule in the rule-based symbology, or as filter on the entire feature.  In both cases, nothing is ever displayed.

So, is what I'm trying to do possible (show/hide features in one layer based on selections in another) and, if so, is my approach flawed or just buggy?


texts: New Error [123] on previously running Python code

texts:
I'm running a toolbox that had previously been able to run smoothly a week prior.  The past two days I've been getting the following error:

Traceback (most recent call last):
  File "<string>", line 301, in execute
  File "C:\Program Files\RCAT-1.0.11 -EPP\VBET.py", line 73, in main
    os.mkdir(projPath + os.path.dirname(DEM) + "/Flow")
WindowsError: [Error 123] The filename, directory name, or volume label syntax is incorrect: u'C:\\Users\\Sabrina\\Documents\\VerdeC:\\Users\\Sabrina\\Documents\\Verde/Flow'

Failed to execute (VBETtool).

Here is the Python code referring to in the code:

 if FlowAcc == None:
        arcpy.AddMessage("Calculating drainage area")
        calc_drain_area(DEM)
    elif os.path.exists(os.path.dirname(DEM) + "/Flow"):
        pass
    else:
        os.mkdir(projPath + os.path.dirname(DEM) + "/Flow")
        arcpy.CopyRaster_management(FlowAcc, os.path.dirname(DEM) + "/Flow/" + os.path.basename(FlowAcc))

In researching this error, most of the time this pops up due to an invalid symbol or character so I've tried changing the project folder to one without potentially problematic symbols or characters (Note this run involved a folder with only one word folder names), but continue to get the same error.

I did recently need to run a repair of my ArcGIS due to issues I was having getting it to run while working on this project and am not sure if maybe that may have affected it.  Previously it was running really slowly, certain licensed tools weren't working and others were not running properly at all. So I'm not sure if that's related or what.  I'm currently working in ArcGIS 10.5 and Python 2.7.1.


texts: Converting distance in ArcGIS Android SDK?

texts:
I have a Polyline that I draw on the map and I would like to get its distance in meters.

right now it's not working, basically what I have is

   distanceSeg=new Polyline();
   double  len1=distanceSeg.calculateLength2D();
   double m_len=
         Unit.convertUnits(len1,Unit.create(AngularUnit.Code.DEGREE),Unit.create(LinearUnit.Code.METER));

The call to Unit.convertUnits throws a GeometryException with "invalid_call".

Any ideas?


texts: ArcGIS Kernel Density with polyline, search radius / bandwidth calculation

texts:
I'm conducting a density analysis of a series of polylines and I'm hung up on determining an appropriate search radius.

The default calculation performed by arcgis is to take the shorter dimension of the input feature extent divided by 30. However, I would prefer to use the formula set forth by Bailey and Gatrell in Interactive Spatial Data Analysis (1995).

They suggest that the bandwidth for a kernel analysis be set to 0.68n^(-0.2) where n=the number of features being analyzed.

This calculation is straightforward for points as each point is the same size. However, the polyline dataset I am working with has lines of varying length. Simply taking the number of lines as n doesn't seem appropriate.

I have considered converting the polylines to points but before I begin with this approach I am curious if there might be other valid approaches?


texts: Displaying representation symbology description in legend of ArcGIS for Desktop?

texts:
When a layer with Representations is turned on, how do you change the 'Description' in the Legend. The setting is normally under the symbology tab and it is not visible when Represenations are turned on. The layer is shown in the legend as the field it is shown on.

Can you modify the 'description' under symbology to show the title of the layer as any text you type in?

I'm using ArcGIS 10.2 for Desktop.


texts: clipping layer using query builder in QGIS

texts: I have a global vector layer which includes several fields in the attribute table. I want to filter one field above a certain threshold and then clip the corresponding areas to make a new layer. I have managed to filter them using the query builder, but I can't find a way to extract those polygons to make a new layer. I'm using QGIS 3.2.