texts:
I am currently creating a vector layer like following:
outLayer = QgsVectorLayer('Linestring?crs='+ crs, 'graph_with_azimuth_ogr' , 'memory')

The outLayer is now stored in memory, which means it only stores the layer temporarily. If I restart QGIS, the layer will be deleted again.
My question is: How is it possible to create a vector layer which is permanently stored or which storage provider do I need to achieve this?
I have already tried it with the "ogr" storage provider, but it does not work. Any suggestions?


texts: Elevation lines label positions (select manually?) in QGIS

texts:
I have a layer of elevation lines in QGIS. Now I want to create a map that shows a certain extent of this layer and shows the labels for certain lines (ideally once for every line).
So far, this isn't a problem. However, I also want to show an inset above the main map. Now I have two problems:

Some labels are printed behind the inset. I can set the frequency of labels along the line from "No Repeat" (label just once per map extent) to a certain distance, however, this will result additional labels that I do not want.
Most lines go beyond the extent of the map I create and come back somewhere else. On the map, this will result in two lines of which QGIS will label only one. However, I cannot assume someone who looks on the map knows that those two lines belong to the same circular object.

Hence, it would be great if there was an option to set the positions of the labels along these elevation lines manually or if there was a trick to get the labels at the right positions.


texts: Python script for an automatic map composer have a bug with refresh

texts:
I have an issue with my script. I would like to show a part of map centered on selected features with a template. The script open the print composer but it don't refresh correctly the map, it's blank. I followed this tutorial. In the map element of the composer, scale is set to 0.0000000000, min X max X min Y max Y are set to 0.000.
here the concerned part of my script:

                self.iface.mapCanvas().zoomToSelected()
                self.iface.mapCanvas().zoomScale(1000)

                canvas = QgsMapCanvas()
                myFile = os.path.join('C:\Users\mypath\\', 'template.qpt')
                myTemplateFile = file(myFile, 'rt')
                myTemplateContent = myTemplateFile.read()
                myTemplateFile.close()
                myDocument = QDomDocument()
                myDocument.setContent(myTemplateContent, False)
                newcomp = iface.createNewComposer()
                newcomp.composition().loadFromTemplate(myDocument)

                composition = newcomp.composition()
                map_item = composition.getComposerItemById('map')
                map_item.setMapCanvas(canvas)

                map_item.zoomToExtent(canvas.extent())

                composition.refreshItems()

My Qgis version is 2.16.3. And the Python console of Qgis don't show any error.


texts:
I am looking for a linux program that allows me to visualize CityGML models (3D). I have tried a whole bunch so far (most of them via wine but also native ones) but I couldn't get a single one of them running.

I am using Ubuntu 13.04 and would appreciate if somebody could answer from personal experience instead of just posting a bunch of links.


texts:
I read the: "Subscript labels in QGIS legends" thread.
I would like to follow up on this for nautical notations.

One Z-value, when displayed, must be displayed in normal script and subscript.
For example: If you want to display the value 9,3 it is supposed to be displayed in the nautical notation as: normal script: 9 subscript: 3

Is there a way to format the label of the value in this manner?
Note that heights above datum have an underscore under the whole digits.

]2


texts: Calculating point distance in separate X/Y-values (lengths of opposite side and adjacent side) using ArcGIS Desktop?

texts:
Description:

I have to position a temporary initial zero point into a corner of my map. Afterwards I want to calculate the geometric distance between many different points and this so called initial zero point.

The challenge is:

I don't need the nearest distance (shortest way) between those points, but the X and Y-components of those distances separately (= lengths of opposite side and adjacent side).

How can I solve that? The nearest distance won't work in my case, because it already shows the hypotenuse distance instead of opposite side and adjacent side lengths.

One possible approach:

Creating a small coordinate system with two lines (X and Y), which are heading through the previous initial zero point. Afterwards I could calculate the distance from those already mentioned points to these lines. As a result, I will get distance-data for X and Y components separately.

Are there any other, maybe more comfortable ways how to do?


texts: How can I combine two datasets with the same schema in to one table?

texts: I have two tables (in CartoDB/PostGIS) that are from two tiger shapefiles that have the exact same attributes, just in different counties. I want to combine them in to a single table, just appending the data from one to the other. What is the SQL call to do that?


texts: Exclude empty variables from OpenLayers 3 popup window

texts:
I have a map displaying multiple locations loaded from a JSON file. The locations are contained in a MySQL database along with various variables for each location including whether or not the location has associated social media pages. Since not all locations have all of the social media pages some of these fields are empty.

Clicking on a location brings up a popup window which displays the location's name with clickable social media icons below it.

I want to exclude the social media icons for empty fields and have tried to write "if" statements to do this but as a Javascript beginner so far have not managed to do this successfully.

I've uploaded my code to JSFiddle but can't get it to display properly there, although my map is displaying correctly on MAMP on my computer (I did try changing the JSON and CSS files to online versions but the map still didn't display properly so I reverted to the original code). The current code shows all the social media icons regardless of whether or not there are entries for all of them associated with that location. This is the JSFiddle link:

https://jsfiddle.net/u5ue682h/5/

Is it possible to place "if" statements within the following "if" statement to filter out the empty variables from the popup?

          if (cfeatures.length == 1) {
        popup_content.innerHTML = '' + cfeatures[0].get('name') + '

 ';
      }


texts: Issue on Removing Index From Field by Arcpy

texts:
Using ArcGIS 10.2 Desktop, I am trying to remove an attribute index from  my shapefile (marten.shp) as:

import arcpy
from arcpy import env

env.workspace = "D:\\Map\\Data"
arcpy.RemoveIndex_management("marten.shp", ["MamerMN"])

Here is the marten.shp Indexes

but I am getting this error 000192

Can you please let me know what I am doing wrong?


texts: Poor Performance when using Spatial Indexes in MySQL

texts:
Re-post of a question asked on Stack Overflow when it was suggested this would be a better forum.

I'm trying a little experiment at pushing a data set which is not geo-spatial but fits it quite well and am finding the results somewhat unsettling. The data set is genomic data e.g. the Human Genome where we have a region of DNA where elements like genes occupy specific start and stop coordinates (our X axis). We have multiple regions of DNA (chromosomes) which  occupy the Y axis. The goal is to bring back all the items which intersect two X coordinates along a single Y coordinate e.g. LineString(START 1, END 2).

The theory seemed sound so I pushed it into an existing MySQL based genome project and came up with a table structure like:

CREATE TABLE `spatial_feature` (
  `spatial_feature_id` int(10) unsigned NOT NULL AUTO_INCREMENT,
  `external_id` int(10) unsigned NOT NULL,
  `external_type` int(3) unsigned NOT NULL,
  `location` geometry NOT NULL,
  PRIMARY KEY (`spatial_feature_id`),
  SPATIAL KEY `sf_location_idx` (`location`)
) ENGINE=MyISAM;

external_id represents the identifier of the entity we have encoded into this table & external_type encodes the source of this. Everything looked good and I pushed in some preliminary data (30,000 rows) which seemed to work well. When this increased past the 3 million row mark MySQL refused to use the spatial index and was slower when it was forced to use it (40 seconds vs. 5 seconds using a full table scan). When more data was added the index started to be used but the performance penalty persisted. Forcing the index off brought the query down to 8 seconds. The query I'm using looks like:

select count(*)
from spatial_feature
where MBRIntersects(GeomFromText('LineString(7420023 1, 7420023 1)'), location);

The data going into this is be very dense along the Y dimensions (think of it like you've  recorded the position of every building, telephone box, post box and pigeon on a very long road). I've done tests of how R-Indexes behave with this data in Java as well as others in the field have applied them to flat-file formats with success. However no one has applied them to databases AFAIK which is the goal of this test.

Has anyone out there seen a similar behaviour when adding large quantities of data to a spatial model which is not very disparate along a particular axis? The problem persists if I reverse the coordinate usage. I'm running the following setup if that's a cause

MacOS 10.6.6
MySQL 5.1.46


texts: Extract sum of 9 pixels into point feature class

texts:
I have a task to create a grid (point feature class)(using Create Fishnet) in the extent of population raster (value of 1px represents number of people living in that area)  and then extract values of the raster into this grid. However, in each point of the grid, I don´t want only the value of the pixel right below the grid point, but sum of values of this point and 8 neighboring pixels, so sum of 9 pixels. I tried using Extract Values to Points, but it only allows to extract the 1 value of the pixel below the point. The idea is that in the end, summing up these extracted values from the points should be equal to summing up values of all pixels of the underlying raster. So any idea, how should I go about this? Is this possible in arcpy?

[The point feature represented by black dot shoud contain attribute value of sum of raster pixel values of 9 neighboring pixels arround and including the pixel under the point (9 pixels that should be summed up are represented by distinct colors]


texts: What is the best way to take a large mosaiced raster dataset (e.g. from http://viewer.nationalmap.gov/viewer/) and tile it to USGS 3.75 minute quarter quad specs?  The problem I am having is creating 300m overlapping borders and applying USGS naming conventions to the output.


texts: QGIS not available in Pop! Linux repo so how do I install it

texts:
It seems the syntax of Pop! OS is such that I can't search this site or google without it being confused as pop up.
It is a scientific based development friendly OS based on Ubuntu.
However you can't get QGIS through a repo so I wondered if anyone knows what the problem is or if there is a work around.
It seems the perfect environment for creating a python extension for QGIS so that is what  I'm trying for.


texts: How to overlay two shapefiles with different coordinate system?

texts:
It seemed to be easy but make big problems, perhaps it is a beginner mistake. Anyway I have to ask you: I need to lay two shapefiles over each other to make a spatial query. The coordinate systems of both shapes are unknown. But I now the one shape was projected in ETRS-1989_LAEA originally; it is a download from here (Berlin). The other shape file was projected in DHDN_Soldner_Berlin.

Now I tried to project the first shapefile to its original coordinate system with the tool “define project” and did the same with the second one. Loading both shape files with different coordinate systems to a data frame the shapefiles are not lying over each other. So I tried to import both to a feature dataset in the File geodatabase. The feature dataset got the coordinate system from the first shape so ETRS-1989_LAEA. Then I import both shapefiles to the feature dataset in the File geodatabase. But this is also not working; the shapes are not laying over each other.

But did I wrong, how can I get both shapes which having different coordinate systems over each other? On the fly is also not working, when I project the coordinate systems from one shape.

Projecting the first shape in the coordinate system from the second shape with the tool “project” is working. But load both shapes having DHDN_Soldner_Berlin as coordinate system to ArcMap 10.2 is again not working.

I read in the forum about the issue but I cannot find my mistake.

Here is the extent from the first shape

Here is the extent from the second shape


texts: I importing live data to CartoDB using Synced table feature which is updated hourly. I also applied simple SQL to the live synced table. However, the selection coming from SQL query would not be updated hourly. Is there any way that I can set the SQL in CartoDB to be updated as well every hour?


texts: Python Hyperlinking Script in Arcmap 10.2

texts:
I have the below python script which recursively walks through a directory and finds, then opens the associated parcel pdf file. How would I modify this to not be recursive but to just look in the directory of "\corp.qr.com.au\app\QRGIS\QRData\Leases\Lease_Plans\Active" I can't seem to get it to work. Thanks

import webbrowser, os
def OpenLink([File_Name]):
  for r,d,f in os.walk(r"\\corp.qr.com.au\app\QRGIS\QRData\Leases\Lease_Plans\"):
    for image in f:
      if image.startswith(str([File_Name])):
        webbrowser.open(os.path.join(r,image))
  return


texts: Are database client files platform specific?

texts:
I'm just looking for some clarification on a question that came up with me this morning.  Are the DLL client files (libeay32.dll, libiconv-2.dll, libintl-8.dll, libpq.dll, and ssleay32.dll) platform specific?

ie - If I have QGIS and ArcDesktop installed on the same machine, can I take a copy of those DLL files from the QGIS installation and move it to the Arc \bin\ folder and then connect to a PostgreSQL database?  OR, do I need to have a specific set of the DLLs for each platform (QGIS/ArcDesktop)?

EDIT - For Clarification:

This ESRI Help Article says that in order to set up a connection between ArcGIS and a PostgreSQL database I need to download a set of DLL files from the ESRI Customer Care Portal and place them in the \bin\ folder of my current ArcGIS installation.

I noticed that I already have a set of these files that were installed with QGIS.  Can I just use that copy (from QGIS) in place of the ones I would download from ESRI?


texts: I have to process 46 MODIS MOD15A2H images (leaf area index) in QGIS 2.12 Lyon. I have to clean the images from bad pixels using MODIS Quality Control. Is there any method to automate the process for this amount of images?


texts:
In ubuntu 16.04 I cannot import gdal in Python. My code starts with:

from osgeos import gdal

but I got error:

    File "ccomt.py", line 1, in <module>
from osgeo import gdal
File "/home/neptun/miniconda2/lib/python2.7/site-packages/osgeo /__init__.py", line 21, in <module>
_gdal = swig_import_helper()
File "/home/neptun/miniconda2/lib/python2.7/site-packages/osgeo/__init__.py", line 17, in swig_import_helper
_mod = imp.load_module('_gdal', fp, pathname, description)
ImportError: /home/neptun/miniconda2/lib/python2.7/site-packages/osgeo/_gdal.so: undefined symbol: VSISupportsSparseFiles

Does anyone faced the same problem?


texts: How to limit Create Fishnet extent by mask & carry over attributes?

texts:
I have a 5km polygon grid of an area which is not a standard square/rectangle, it is in fact a rough ovoid made up of around 500 5km2 cells.

What I want to do, eventually, is have a 100m polygon grid covering the entire area that is currently represented by the ovoid(ish) 5km grid - about 1.25 million cells.

The caveat is that the 100m polygons should be linked to a 'parent' 1km grid square: i.e. you have 12,500 1km cells with a given ID (1KM_ID) and then the 100m cells are given an ID too (100M_ID) which goes from 1 to 100 for each 1km cell.

Using Create Fishnet will create 1km cells based on the bounding box of the 5km cell ovoid, which will greatly increase processing time, i just want the fishnet based on the cells that exist. Same applies for then creating 100m cells from the 1km cells.

This can be done really nicely in FME but for some reason there are processing issues and duplicates are created around the edge of each 1km cell, messing the whole thing up. Any ideas for AcrMap?


texts: This seems like a simple question, but I'm having trouble finding out the answer. I'm building a map in tilemill and using data from a geojson file. I plan on hosting it later using tilestream. Once I do, will I be able to update the underlying data and have it be reflected in the map? Or will I need to reexport and deploy again?


texts: How to close all tooltips defined against the map object?

texts:
In Leaflet it's possible to bind one or more tooltips to an individual marker, or a layer, or the base map object.

How can I close all open tooltips that have been bound to the map object?

var map = L.map("map");
L.tileLayer("http://{s}.tile.osm.org/{z}/{x}/{y}.png").addTo(map);
map.setView([48.85, 2.35], 12);

// Define some tooltips against the map object
var tooltip1 = map.openTooltip("tooltip 1",[48.84, 2.34]);
map.openTooltip("tooltip 2",[48.83, 2.33]);
map.openTooltip("tooltip 3",[48.82, 2.32]);

// Attempt to close tooltips, individually or en masse - neither works
map.closeTooltip(tooltip1);
map.closeTooltip();

See also the JS Fiddle version of the above code.


texts: How to use in_field_map parameter of CreateCompositeAddressLocator_geocoding?

texts: I would like to know if anyone who has worked with the CreateCompositeAddressLocator_geocoding method knows how the "in_field_map" parameter is formed?


texts: Walking a number of kilometers in one direction in digital map

texts:
Is there a recommended way to calculate the new position of a point in a map when "moving" the point for example 10 kilometers to south ?

Example:

Point = {
    lat : 'LATITUDE VALUE',
    lng : 'LONGITUDE VALUE'
}

Now we want to get the new Position of the Point when the Point is moving like this:

Point_2 = movePoint(Point, 10, 'kilometers', 'south'); // returns lat,lng

Now I am looking for a recommended way to calculate something like this or maybe with a framework like turfjs (I think turfjs is not providing this ?)


texts: Why am I unable to publish a postgis layer to geoserver?

texts:
I left arcgis and microsoft sql and now I am trying to connect Qgis+PostgreSQL (PostGIS) + GeoServer.

I've brought my postgis database to the server, created a "store". And now I'm trying to publish a layer. So I went to Layers- Choose my layer- Publish- Edit Layer (Data/Publishing)- configure everything there- SAVE and nothing happens... I can't see my layer in "Layer Preview"

So maybe configuring was wrong. I chose the native SRS, point style, authority name: admin, authorithy URL (is it the server's URL?? ). What is layer identifier, as well?

Thank your patience and answers:)


texts: What to do with gigantic gdb directory, looking for snaller geojsons

texts:
I was given a large .gdb directory containing all parcels for an entire state. I tried converting the entire thing using ogr2ogr using ogr2ogr -f GeoJSON TEST.json FOO.gdb, but the resulting GeoJSON file was 35.5GB, then the process died with Terminating translation prematurely after failed.

Is there a notion of decomposing a .gdb file into smaller files, so that I may try to convert them into parcels in geojson format, say per parcel, or per town maybe?


texts: Remove all labels from Esri basemap for use in server run application

texts:
We pull in basemaps into our mapping tools, but we would like to remove all labels, not just the reference layer's labels of a basemap. I'm struggling with how to achieve this without having to build a basemap through reverse engineering. This basemap:
https://www.arcgis.com/home/item.html?id=a2eae30d60e34cafaff0cc28f8edc324
is the closest I've found, but I can't see a way to add this item to my ArcGIS Portal/Enterprise so that I can add it to our tools. What am I missing?


texts: Adding editable Excel (or Access) table in ArcMap?

texts:
I have flood data from the past ten years written in an Excel table. Now I would like to connect the table to ArcMap, so I can join it with the shapefiles. Normally this isn't that hard, but because we talk about floods, we want to hold the Excel table up to date.

What we want is a connection between ArcMap and Excel/Access so we can edit the tables and the data will automatically change in ArcMap (without a new import->join->etc..).

I think it is something with a geodatabase, but till now it isn't working at all.

Does somebody know if this is possible and how?


texts: How to create a loop in R to average ndvi raster time series?

texts:
I have monthly ndvi raster geotiffs for multiple years. I need to run a loop to get the average of each of these ndvi grids/rasters.

  Each of my ndvi rasters contains 143 pixels. I would like the mean of all the pixels within a raster across my time series.

Can someone offer any advice?


texts:
How can I get spatial statistics on my vector (tab) data in MI 15 64-bit. I just need to see total line lengths, areas etc as appropriate.

From the MI help I tried the following

  To see the details of a data record:

  Open a .TAB file with statistical data in it. Select a geographical
  object in the Map window. On the HOME tab, in the Tool Windows group,
  click Statistics.

But all I get is

  Table CCRC_Road_Centerline_NEW has 464 records selected.


texts: Which coordinate system will prevent eastern Asia from splitting?

texts:
I would like to create a simple graphic for a publication showing all continents (without Antarctica) but I would like Alaska to stand alone and the very eastern tip of Asia to be attached to the main Asia.
I like the way Sphere_Gall_Stereographic WKID: 53016 Authority: ESRI looks but like WGS84 GCS this CS like many others splits the tip of Asia from the main continent.

In some maps (also in some North American news shows on TV where I noticed this), Asia in fact appears as a single continent and Alaska stands alone.
Although this may be technically incorrect, visually, I like the continuity.

I tried Robinson projection (both world and sphere) which is what is used in this World Map but whereas the linked map shows Asia continuous Robinson in ArcGIS splits it apart. I am confused.

Is this something I can Achieve in GIS by using a specific or modified coordinate system or do I just have to export the map to a vector graphic and attach these two pieces manually in post?

Sphere_Gall_Stereographic ESRI WKID: 53016

WGS84:

Robinson in ArcGIS:

Robinson (World Political Map, Wikipedia):


texts: How to solve "error 000816: The tool is not valid" with ET Geowizards tools?

texts:
I have installed ET Geowizards 11.2 for ArcGIS 10.3 using setup.exe and followed the instructions. I am using ArcGIS 10.3.

I can use ET Geowizards tools from the toolbar and I can run them in ModelBuilder. However, I cannot run them from the toolbox. When I load the ET Geowizard.tbx into the ArcToolbox and run a tool, I get the following error:

According to online resource http://www.ian-ko.com/ET_GeoWizards/UserGuide/ToolBox/toolbox_howto_load.htm the tools get registred in the installation. However, this doesn't appear to have been succesful.

My questions is how to solve this problem.


texts: Google Maps zoom number and min/max scale(Openlayers)

texts:
How to set Google maps' zoom number accordingn to min/max scale of a wms layer? If min scale is 1/50000 and max scale is 1/5000 (units: m), which zoom number is resonable? I tried to set min zoom number to 10, it worked fine, but I got problem with max zoom number.I tried 20 and 22, but Google Streets and Google Physical gave no responce when zoomed in (only the wms layer zoomed, see the screenshot).

Any tips how to set zoom number?


texts: Docker Image: Atmospheric Correction in Earth Engine

texts:
Following https://github.com/gee-community/ee-jupyter-contrib/tree/master/docker/atmcorr-ee I was able to build and run docker. And to authenticate and git clone, but while:

cd gee-atmcorr-S2/jupyer_notebooks/
jupyter-notebook sentinel2_atmospheric_correction.ipynb --ip='*' --port=8888 --allow-root

then I got:

root@9c6555522c26:/gee-atmcorr-S2/jupyer_notebooks# jupyter-notebook sentinel2_atmospheric_correction.ipynb --ip='*' --port=8888 --allow-root
[I 20:18:29.195 NotebookApp] Writing notebook server cookie secret to /root/.local/share/jupyter/runtime/notebook_cookie_secret
Traceback (most recent call last):
  File "/miniconda/lib/python3.6/site-packages/traitlets/traitlets.py", line 528, in get
    value = obj._trait_values[self.name]
KeyError: 'allow_remote_access'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/miniconda/lib/python3.6/site-packages/notebook/notebookapp.py", line 864, in _default_allow_remote
    addr = ipaddress.ip_address(self.ip)
  File "/miniconda/lib/python3.6/ipaddress.py", line 54, in ip_address
    address)
ValueError: '' does not appear to be an IPv4 or IPv6 address

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/miniconda/bin/jupyter-notebook", line 11, in <module>
    sys.exit(main())
  File "/miniconda/lib/python3.6/site-packages/jupyter_core/application.py", line 266, in launch_instance
    return super(JupyterApp, cls).launch_instance(argv=argv, **kwargs)
  File "/miniconda/lib/python3.6/site-packages/traitlets/config/application.py", line 657, in launch_instance
    app.initialize(argv)
  File "</miniconda/lib/python3.6/site-packages/decorator.py:decorator-gen-7>", line 2, in initialize
  File "/miniconda/lib/python3.6/site-packages/traitlets/config/application.py", line 87, in catch_config_error
    return method(app, *args, **kwargs)
  File "/miniconda/lib/python3.6/site-packages/notebook/notebookapp.py", line 1628, in initialize
    self.init_webapp()
  File "/miniconda/lib/python3.6/site-packages/notebook/notebookapp.py", line 1378, in init_webapp
    self.jinja_environment_options,
  File "/miniconda/lib/python3.6/site-packages/notebook/notebookapp.py", line 159, in __init__
    default_url, settings_overrides, jinja_env_options)
  File "/miniconda/lib/python3.6/site-packages/notebook/notebookapp.py", line 252, in init_settings
    allow_remote_access=jupyter_app.allow_remote_access,
  File "/miniconda/lib/python3.6/site-packages/traitlets/traitlets.py", line 556, in __get__
    return self.get(obj, cls)
  File "/miniconda/lib/python3.6/site-packages/traitlets/traitlets.py", line 535, in get
    value = self._validate(obj, dynamic_default())
  File "/miniconda/lib/python3.6/site-packages/notebook/notebookapp.py", line 867, in _default_allow_remote
    for info in socket.getaddrinfo(self.ip, self.port, 0, socket.SOCK_STREAM):
  File "/miniconda/lib/python3.6/socket.py", line 745, in getaddrinfo
    for res in _socket.getaddrinfo(host, port, family, type, proto, flags):
socket.gaierror: [Errno -3] Temporary failure in name resolution

Any hint on how to proceed?


texts:
I am totally new at this field and I want to use Rinex data to add more precision to my GPS receiver.

I don't know where to start and even I spend a lot of time googling, I just got more confused.

Can anyone help me with a starting point?

I want to determine the exact 3D position of a flying quadrotor while flying (So i need the position realtime).

The main project is multi-agent UAV flying, so their position has to be accurate enough to prevent crashing during flight. Also for accurate landing in a desired position.

In the first phase I need to know position relative to a ground station in a not so far distance (I mean about 5-7 meter from robot). And this station (the positioning reference) doesn't move.


texts: How to render a vector layer using a QgsMapCanvas in QGIS 3?

texts:
I updated this ShapeViewer Demo for QGIS 3 and the layer does not render. All that shows is a white rectangle.

The updated code is here. What am I missing?


texts: ArcGIS Server - How to give different permissions to different users or roles?

texts:
I have a map service with Features Access capabilities enabled on ArcGIS Server 10.3. This map service is protected using ArcGIS Server tokens.

I would like to have two categories of users with different permissions:

Users Category 1: Users that can create/delete/edit/query data using Features Access REST API. For this I use the Editor Tracking
Users Category 2: Users that can only query specific data

For Instance:

User_A_Editor (Category 1)
User_A_Viewer (Category 2)

User_A_Viewer would be able to only query data created by User_A_Editor. User_A_Viewer won't be able to create any data.

I know there is a possibility with Features Access coupled with Editor Tracking to restrict editing/deleting/query to features created by the user, but it seems I can't prevent a user of the Category 2 to create new feature or to query data created by any other users of the Category 1.

I was thinking to create a Geoprocessing Service that will be used as a kind of proxy for any REST queries made to this map service. So in the client application, all queries would be made to the GP service instead. This GP will verify if the user is allowed to do such REST query, and if yes, the GP will call the REST API and return the result to the client. But in reality, that won't prevent a wise user to use the REST API in a way not wanted as it will always be accessible.

Is there any solution?


texts: Re-rendering tiles after updating database using mapcache_seed

texts:
I have developed the map using MapServer and pre-generated tile using mapcache. Now I have an idea to update the database every day. As changes apply in the database only those parts should be re-render and new tile should be generated. As changes are very small and rendering all ties from start is difficult to process.

I got expiry list of tiles by following osm2pgsql:

osm2pgsql --append --slim --cache 1000 --multi-geometry /home/parsis/diff_file_minute.xml --number-processes 2 --hstore -d testdemo -r xml -e6-11 -o /home/parsis/expire.list**

Got list something like this:

 - 11/329/713
 - 8/40/91
 - 10/164/357
 - 7/29/47

I find another approach using render_expired, but don't know how to use it with mapcache.

How is it possible using mapcache_seed?


texts: How to clip buffer polygon from its source line?

texts:
I'm trying to clip a 30 meter buffer polygon built around a highway line feature using a smaller 6 meter buffer that represents the actual width of the road, so that the original 30  meter buffer is actually 6 meters away from the highway line feature it was created from.

This is so I can calculate slope and land cover, etc., on the edge of the road without including the road itself.

Whenever I try this, it just clips the 30 meter buffer down to 6 meters.

Is it possible to clip a buffer off of its source feature in this way?


texts:
I am trying to figure out if there is a way to multiply two different sets of info (in different columns and cells) and have the answer end up in another column/cell in the same feature. But not only multiply those values for data references, but to become a functional formula so that if the two sets of values change, the answer will change automatically? And if this is possible, can it be shared to arcgis online.

I would like to multiply acres and manure application rate and it be answered in the Total Manure Column, like in the row "Adams Farm". I know how to calculate it after the data is entered. But I would like it to be automatically answered if the either acres or manure application rate changes.

I would also like this to be uploaded to arcgis online and used with the collector app. Is this possible?


texts: Combining spatially identical polylines in the same layer

texts:
I have a shapefile of road volumes as polylines, and want to be able to show the sum of the volumes for polylines that are spatially the same/similar but are for travel in opposite directions. For e.g. one polyline for traffic travelling North may have a volume  of 8000, while the equivalent polyline that is for traffic going South has a volume of 6000. I need to be able to locate each of these polylines and sum their volume values, and show this as a separate layer. Note that some polylines within the dataset are along one-way streets and therefore don't have an equivalent polyline that travels in the opposite direction.

I've been able to distinguish each line pair in the attribute table by generating the xmin and xmax, however I'm not keen on selecting each pair manually and calculating the sum of volumes given there are a few thousand. I've tried spatial join (MMQGIS plugin) and union in QGIS but haven't had any luck. I've also toyed with the QGIS Field Calculator a bit but I'm not very familiar with how to use some of the geometry commands for this type of analysis.


texts: How to find disconnected islands in a road network layer, using QGIS

texts:
I have a road network layer that can't route from one side of the layer to the other.  Topologically, there is no connectivity between the start and end nodes, probably due to missing links.  This creates disconnected "islands" in the network graph.

How can I use QGIS or GRASS or a python plugin to analyse my road network, perhaps traversing it from the starting node, to the extents of the connected links?

The output can either be a group column, where all connected links will all have the same ID, and the number of distinct IDs will be the number of disconnected islands, or a colour coded equivalent, like the one below, produced by TransCAD.

I am using QGIS 2.14.0-Essen.  I have tried Topology Checker, Processing modules (formerly Sextante toolbox), GRASS v.net, v.clean, etc, and PyQGIS Developer Cookbook (http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/network_analysis.html#areas-of-availability)

The next step is to try and iterate, from the starting node, all "touching" links, adding them to a selection set, until there are no further ones to add.  Label these as group 1 and then move onto another link that has not been selected and repeat for group 2, and so on, until there are no more unselected links.  But how do I do this easily?

Other tools, like ArcGIS, have a "Find Disconnected" tool.  But I am looking for a QGIS solution, please.  An alternative might be the pgrouting.org extension for PostGIS.


texts: Passing output stored in single variable into another script using ArcPy?

texts:
I am running set of processes each with different outputs in a python script.

Is there a way to set a variable as global so that just the output of that single step can be accessed from other scripts?


texts: ArcGIS geocoding suggestion service, max suggestion number

texts:
How do I set the maximum number of suggestions returned by ArcGIS suggestion service to more than the default 5?
In the docs (docs) they say,

  "maximum number allowed by the service".

So I believe, that this is a thing you could change.

I am using ArcGIS Server 10.3.1.
Here is a sample from the docs with more than 5 suggestions link


texts: Choosing CRS to use for South Africa that will give map scale in meters in QGIS?

texts:
This has been a constant source of frustration for me. I can't seem to always get my map scale in meters in QGIS, even if I use a CRS that has map units of meters, e.g. Alberts Equal Area and others. Sometimes it works, but most of the time it doesn't.

I usually work with an area of a few square kilometers. I then have to georeference a satellite image of the area and load my different vector and raster data.
Let's specifically look at the georeferencing, I first deselect OTF, then I select my GCP points and set my transformation type. Then I select the CRS, e.g. Hartebeeshoek94/Lo25 (map units = meters) and run the georeferencer. The result is a georeferenced satellite image with the correct coordinates but a scale of 1:1. The scale bar is then also completely incorrect in the map composer.
I have now found that every time I save the georeferenced image as a CRS with the units in meters, the scale loads as a 1:1. When I load the layer as a CRS with map units in degrees, the scale is correct by the scale bar is still incorrect in the map composer (as would be expected).


texts: Ratio of area between park and residential parcels

texts:
I have three shapefiles: park, residential parcel, and participants of a survey (single Point shapefile with 100 dots).

I want to find a ratio of area (park/residential parcel) within 1 and 0.5 mile buffer of the applicants. I have total applicants of about 100.

What are the steps that I need to follow to perform this work in ArcMap?

So far I created 1 and 0.5 mile buffer around the applicants and clipped the park area and parcel areas with the buffer shapefile.

What should I do next so that I can find the ratio of area used by each applicants?

I am using ArcMap 10.2.2. Is it possible through ArcMap?

In the figure, two red dots are my applicants, the blue colored polygons are clipped parks within 1 mile buffer of those two applicants, and light apple green color is parcels clipped within those buffer. Now i need a ratio of area between those blue color to light apple green color within each buffer.


texts: Is it possible to produce a decent looking shaded-relief map at 50K-scale-or-larger using SRTM3 DEM?

texts:
I'm trying to make a decent-looking shaded-relief map of a particular area in 50K-scale-or-larger using open-source tools. The best (and free) option I have for DEM data for that area right now is the global 90m resolution SRTM3 data.

So far I have tried using GDAL tools (gdalbuildvrt,gdalwarp,gdal_translate,gdaldem) and GRASS tools (r.resamp.rst,r.shaded.relief) to produce the shaded-relief. But the output always retained these characteristics:

Image is too pixelated at 50K-scale-or-larger. I've tried re-sampling them to 30m and even 5m using RST/bilinear/cubic methods, but the images are somehow still pixelated
The landform is not clear or too generalized. I have to increase the z-factor up-to 4.0 to obtain better landform. I don't know if this is normal..

I've read most of the work done by Tom Patterson in http://www.reliefshading.com and http://www.shadedrelief.com, and I'm trying to emulate his work - but I realized that most of the pointers he gave revolved around small-scale maps, and the use of propriety softwares to get the job done (I can only afford to use open-source stuff for now). So I'm at a lost right now.. Perhaps the enhancement for the SRTM3 DEM can be done in GIMP somehow?

So my ultimate question: is it possible to make use of SRTM3 data to make a decent 50K-scale-or-larger shaded-relief maps? If I can, how should I go about processing it?

Appreciate any leads on this matter, thanks!


texts: Converting ascii windspeed files to rasters

texts:
I have been trying to convert windspeed files to rasters without success. It was discussed in an Arcgis forum which I followed. This involved replacing semi-colons with spaces and adding a header.

I was hoping to include one of the original .asc files together with a .csv version I 'cleaned up'.  Unfortunately I cannot see how to do this without copying and pasting the text.
The first line of the original .asc file is:

Conversion of speed45i.dat   Wind m/s  Wind Speed @ 45m N. Ire
(100,  499); 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0;

I changed this to:

ncols 700
nrows 1300
x||corner 0
y||corner 0
cellsize 1000                                                                                "
"( 0 1299) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;"

When I usen the raster translate convert format plug-in I got the error message:

  The process failed to start. Either the invoked program is missing, or
  you may have insufficient permissions to invoke the program.

Any suggestions gratefully received.


texts: Estimating the water coverage area behind a dam

texts:
I have a task to estimate the area that will be covered by water when a dam is built in the provided location. The only thing I have is the max water level calculated for the reservoir.
Is that enough or I will need more details and data?
Can it be done just using GIS tools?


texts: Copying graduated scale from polygon to point layer in QGIS?

texts:
I have a polygon layer (represented by cells) and a point layer. I've used the graduated symbology with the polygon layer.

Now I'd like to export this graduated scale to the point layer. The field of classification has the same name in both layer.

I've tried to copy and paste the style but it didn't work (probably because I'm trying to export something that represent a polygon to a point).

Am I missing something?


texts: How to Apply a Mathematical Equation on a Feature Collection in Google Earth Engine?

texts:
I have uploaded a table in GEE with height and diameter columns.Now I want to write a function that would apply an equation to each feature. I am unable to use a simple equation as I might have done in excel. The '.expression ' in GEE works only on images not on feature collection. Please suggest a simple method to do it.

Here is the link to my fusion table: https://www.google.com/fusiontables/DataSource?docid=1tuVLITArmX663I8wGZ3l6h-hAN9_q6wO_5papUwi

This is the equation I want to implement :  ln(AGB)=5.349+0.96*ln(D^2H)

where D is for DBH (m) and H is for height (m)


texts: Why does the transformation bring bad results?

texts:
I have this map which has a polygon drawn on it. For the background map I use a WMS via ol.source.TileWMS and the polygon uses a ol.source.Vector object. Normally my view projection is EPSG:31467 but now I tried to load another WMS which does not know  EPSG:31456 but EPSG:25832 instead.

I know up to day OL3 cannot reproject raster data, so I decided to automatically create a new view-obj when the new WMS (more precisely the layer) becomes visible. The new view has the transformed (from 31467 to 26832) center of the old view. As the polygon is in the wrong place now, I transform its geometry as well. Eventually I experienced the problem that the view’s center and the polygon seem to be translated by about 40m.

Why is that? I do not care for the visual difference between the different center positions, but the polygon is in the wrong position now. It is in the same position on screen but not on the map. I thought it should work when I create a new view, and I do not expect the transform function to make such a big difference.


texts: Joining records from two tables that are not completely identical

texts:
I am trying to join a table and a shapefile using ArcGIS 10.2.1 and I need some help. The table has simple records similar to this one:
OID1 WORK1
1   XW-02600
2   LL-02608
3   PM-02609
4   GL-02610
5   GL-01894
6   XW-01895
7   PM-01884

The shapefile has some concatenated records:

OID2 WORK2
1    XW-02607 / GL-01894
2    PM-02609 / GL-01894
3    PM-01896 / XW-01895
4    GL-02610
5    LL-02608
6    XW-02607 / GL-01894
7    GL-02610

The table and shapefile I have is much more complex than the example above, the table has much less records than the shapefile. The shapefile has many records with concatenated values.

All the records in the table should be in the shapefile (but not vice versa), a simple join will not join anything after the "/". I would like for the join to be able work regardless of being behind the "/".

Is this even possible?


texts: Creating one (png, pdf, etc.) map per feature type from QGIS Composer thematic map

texts:
I've got a shapefile that contains thousands of features that can be group by categories.
These features are sometime located in the same place so they overlap.

An example of the content of my shapefile is :

id | scientific_name
1  | foo
2  | foo2
3  | foo2
4  | foo

I would like to show the location of each feature depending on the scientific name.
So, I've made a thematic map with a different colour for each scientific name. Now, I would to be able to have one map (png, pdf, etc.) per scientific_name (some kind of thematic export).

Using the previous example, I will have :

foo.png (with only the features 1 and 4 visible and hide the others)
foo2.png (with only the features 2 and 3 visible and hide the others)

To give a better explanation of what I would like to do, here are some screenshots of what I'm trying to achieve :

1) this is my main QGIS project. My data are classified with a different colour for each scientific name (all the data are inside the same layer, but I can split them as well in different layers if needed)

2) For each scientific name, I would like to export a map that only display one scientific name

Do you have any idea how I could do that ?


texts: Creating 1 hectare Polygons given latitude longitude in Python

texts:
I'm trying to work with openweathermap's Agromonitoring API.

To use their API you have to create polygons representing your area, then you can request weather data, soil moisture etc for that particular area/polygon. There are 2 ways of creating polygons:

Use their graphical tool to draw the polygon on a map manually.
Use their polygon API to create the polygon given a list of coordinates.

I'm interested in the second method, simple get requests with the following body accomplish polygon creation for me.

body = {
"name":"TestPolygon",
"geo_json":{
   "type":"Feature",
   "properties":{

   },
   "geometry":{
      "type":"Polygon",
      "coordinates":[
         [[0, 0], [0, 0.0009], [0.0009, 0.0009], [0.0009, 0], [0, 0]]
      ]
   }
  }
}

Their API requires a minimum polygon area of 1 hectare, so I tried this request with the above dummy coordinates and the polygon created had an area of 1.003 hectares. I want it to be as close to one as possible.

Now the actual question, is how do I create an appropriate coordinates list given a single point (lat, lng), which would result in a polygon of 1 hectare? I only want to input latitude and longitude, and then get a polygon as close to 1 hectare as possible. Perhaps the single point can be the center of our polygon (and we expand around it), or it can be a corner point of a rectangle shape. The shape doesn't matter exactly, it just needs to be nearly 1 hectare.

Since, my coordinates list brought me a 1 hectare polygon with 0.0009, I thought I could fix 0.0009 to be my factor for all coordinates, like so:

# I tried creating a rectangular polygon.
lat = 25.46666
lng = 68.71666
factor = 0.0009
lat_lng_list = [[lng, lat], [lng, lat + factor], [lng + factor, lat + factor], [lng + factor, lat], [lng, lat]]

but to my surprise, this factor does not remain constant for all latitude,longitude points, And it needs to be changed (through trial and error) until I reach a 1 hectare polygon for that specific point. It's also worth noting that the only way I find out if my polygon is 1 hectare or close to it, is through the response I receive to my get request.

GIS is not my field, so I apologize if this question is silly, but I'm interested in hearing how you would solve this problem. Or how I can use code to automate polygon creation (so the user won't have to manually draw it on their map).

Just to reiterate the problem statement: Given a single point (lat, lng), create a 1 hectare polygon of any shape with that point as the center or as the corner (whichever is easier I have no clue).


texts:
I have a shapefile that consists of 375m x 375m cells from the State of California, consisting of a total of 3,022,037 cells. The raster file covers the same geographical area, but consists of smaller cells sized 50m x 50m. I am trying to run raster zonal statistics to infer aggregate features for each polygon in the shapefile.

QGIS fails to complete this on a 2.5 GHz Intel Core i7 with 16 GB 1600 MHz DDR3 RAM. It freezes; the activity monitor shows that it goes to a non-responding state.
A possible source for this problem is mentioned at Zonal Statistics Freezing QGIS 2.2. However, I ensured that the raster does not have any NA values, and there are no "complex" polygons in the shapefile.
I also tried this with rasterstats, which abruptly ended after an hour. I am currently debugging to check what might have caused this.
I also tried to explicitly code up the lookup in python myself, but naturally, this is significantly slower. It takes about 13-14 hours to work on smaller raster file.

My question seeks to know the standard procedure to deal with zonal statistics when the raster and vector files are large. This process needs to be repeated for 40-50 raster files, so this is a problem. A possible solution is mentioned in the paper "Distributed zonal statistics of big raster and vector data", but I am not sure if there is an open source code that can do this. Any guidance on this would be helpful.


texts: Using arcpy to separate surface objects mixed into one file into different shapefiles

texts:
I'm trying to do some elementary tasks with ArcPy and ArcGIS, namely separate some information in a file into separate shapefiles based on values for a field.

My initial file has a field called "Source" and "Source" contains combined information for each record about its state and type of feature (rail, road, utlility etc.), given in a string format like "NYRDLINES" (NY, roads) or "MAUTLINES" (Massachussetts, Utilities). I want to create shapefile outputs with just information for, say, roads in NY or utilities in MA.

To do this, I've used "AddField_management" to add 2 new fields to the table: "State, and Feature type", but don't know how to retrieve the information from "Source" to separate into my 2 new fields.
I think I should be using "CalculateField_management" to calculate the field values, but am stumped when it comes to defining the expression to be used.

How can I do this?


texts: Extracting land cover data (in percent) from shapefile in QGIS?

texts:
I have two shapefiles, one containing stream catchment areas, the other information on land cover classes. The problem: I need to know how many percent of the catchment area is covered by forest, grassland and houses.

I read some tutorials and I understand that I need to clip the land cover shapefile so that it fits the catchment shapefile, then I can somehow calculate the percentage of each land cover class.

However, no matter what I try, it does not work. I get empty shapefiles all the time. Could someone explain how to extract the land cover data? I am using QGIS 2.8 on a Mac.

They are in the same projection system.

This is land cover information

And this the catchment information

I rasterized the land cover shp and used the stream catchment shp to clip it as you suggested. However, the final layer never includes all the relevant land cover information, but only (for example) the "VEG G5" info. I guess I'm doing something wrong when rasterizing the land cover shp or when clipping the layers?


texts: How to change the projection of a UTM-georeferenced raster in Quantum?

texts:
I have a project file in Quantum which contains some shapefiles and rasters that are properly georeferenced. However, I have started using GeoTiffs which are in UTM - WGS84 that don't seem to be projecting properly.

When I import the GeoTiffs, Quantum reads their X,Y (Eastings, Northings) coordinates like they were lat,lon. How can I convert these X,Y values into lat,lon coordinates? I have set all the layers in the project file to the same coordinate system and I have also enabled on the fly CRS transofrmation to no avail.

Thanks for you help!


texts: Export composite image from Google Earth Engine to Drive

texts:
I am new to GEE and have tried to understand my issue by looking through the existing questions related to exporting images to Drive. Unfortunately, I cannot seem to address my issue and hope that someone may be able to help me.

The problem: the file that is exported to Drive is only 1GB, so there is something going wrong (the file should be much bigger!)

My goal: export the composite image for my study area to GEE, so that I can then open it in Idrisi or ArcMap. I believe that the study area is so large (2million square hectares) that when I export it, it should be exported as multiple tiles in geotiff format. however, I do not know what I am doing wrong in my export code because the output in drive is just a single, 1GB file.

The code I am using looks as follows:

// Define a region of interest  (i.e. the shapefile you imported)
var roi = K2C_aoi;

// Load the Landsat 8 scaled radiance image collection.
var landsat2013 = ee.ImageCollection('LANDSAT/LC08/C01/T1')
    .filterDate('2013-05-01', '2013-07-31')
    .filterBounds(roi);
print(landsat2013);

// Make a cloud-free composite.
var composite = ee.Algorithms.Landsat.simpleComposite({
  collection: landsat2013,
  asFloat: true
});

// Clip composite to aoi... 'clipped_composite' is the image you ultimately want to export
var clipped_composite2013 = composite.clip(roi);

//visualize the clipped composite
//Map.addLayer(clipped_composite2013, {bands: ['B4', 'B3', 'B2'], max: 0.5, gamma: 2}, 'Natural color', true);
Map.addLayer(clipped_composite2013, {bands: ['B5', 'B4', 'B3'], max: 0.5, gamma: 2}, 'Color infrared', true);
Map.centerObject(K2C_aoi, 9);

// Select the image and the band combination.
var compositeSelection = ee.Image('clipped_composite2013')
  .select(['B4', 'B3', 'B2']); //natural colors

// Create a geometry representing an export region.
var geometry = K2C_aoi;

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: clipped_composite2013.float(),////////////////////////remove this?
  //image: clipped_composite2013.mosaic(), not a function
  description: 'imageToDriveExample2',
  maxPixels: 3784216672400,/////////////////////////////remove this?
  scale: 30,
  region: K2C_aoi,
  maxPixels: 3784216672400,///////////////how many pixels should I list?
});


texts: Javascript API: queryTask error - where should I put the alert?

texts:
Where in the code do I add an error message if nothing is returned?

     function execute() {
     //
     var queryTask = new QueryTask("http://sepa-app-gis01/arcgis/rest/services/live/SEARCH/MapServer/2");
     var query = new Query();
     query.returnGeometry = true;
     query.outFields = [
       "POSTCODE"
     ];
     query.text = dom.byId("Postcode").value.toUpperCase();
     queryTask.execute(query, showResults);
 };

Should I add Error to the queryTask.execute first??

queryTask.execute(query, showResults, Error)


texts:
I am using gdal_translate -b ... in.sid out.sid to remove excess bands (there are 5 bands) from .sid files as a layer created from these .sid's will be a JPEG layer. When executing command on a .sid file and expecting to get .sid file as an ouput but I get a .tiff and a warning

Warning: The target file has a 'sid' extension, which is normally used by the MrSID driver,
but the requested output driver is GTiff

"The following format drivers are configured and support output" list does not show MrSID driver among supported drivers.

How to make gdal_translate support MrSID driver?


texts: OpenLayers measure and draw controls don't show the sketches on move in touch devices

texts:
I have a GIS API that uses the OpenLayers 2.13 measure and drawFeature controls. But let's just talk about the measure control, because the issues are certainly connected.

As I show in this CodePen, the measure control is supposed to show the immediate measures, i.e., print them on mouse move.

measureControl = new OpenLayers.Control.Measure(OpenLayers.Handler.Path, measureControlOptions);

map.addControl(measureControl);

measureControl.events.register('measure', measureControl, function(e) {
    printOutput(e);
});

measureControl.events.register('measurepartial', measureControl, function(e) {
    printOutput(e);
});

Of course, when I use this example in a touch device, I can't use the onmousemove event. Instead, the behaviour is like this:

1st tap to start measuring
2nd tap to stop measuring

It shows the desired length, but I am not able to perform an accurate measurement because I want to see the measures as I move my finger on the touch screen.

I am not very familiar with touch events, so I apologize if this is a stupid question. but is there an easy way to mimic the mouse behaviour on the touch device?

The desired behaviour is something like:

tap down --> starts measuring
tap move (without lifting the finger) --> keeps measuring and show the new measurement on each move
tap up --> stops the measurement


texts: Saving output of MMQGIS Geocoder to British National Grid receive error message "Latitude or longitude exceeded limits"

texts:
I have used MMQGIS geocoder multiple times and then saved the output with a CRS of British National Grid (BNG) 27700 with no issue. However, I am currently trying to achieve the same result but get the following error message:

Export to vector file failed.
Error: Failed to transform a point while drawing a feature with ID '48'. Writing stopped. (Exception: forward transform of
(-2.026525, 0.783758)
PROJ.4: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 +to +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
Error: latitude or longitude exceeded limits)

The error states that the latitude or longitude exceeded limits.

When I open up the attribute table for the MMQGIS output I cannot see any obvious problem. There is no x or y fields but that is normal for MMQGIS because it does not create them.

The project CRS is set to BNG 27700 and the output from MMQGIS is ESPG:4326, WGS 84.


texts: Batch Geocoding in R over daily limit with pay as you go?

texts:
I have found the the geocode function in ggmap to be very helpful in the past.  Now I am trying to geocode about 60,000 addresses. The google maps API daily limit is 2500.  Google offers a "pay as you go" service that increases daily limits.

https://developers.google.com/maps/documentation/geocoding/usage-limits

Here is my question:

Is there a way to access a "pay as you go" project through the ggmap geocode function (or something similar)?  If so, how would you do that?  That is:  how to use R's "geocode" for batch geocoding where you are willing to pay to go beyond the daily limit?


texts: How to do regression analysis out-of-memory on a set of large rasters in R?

texts:
I am trying to do a regression analysis on a set of large rasters (254,004,000 cells each).  Ultimately, I want to run something like the following (or a bit more complex, but let's start simple!):

model<-lm(dv ~ iv1+iv2+iv3... data=df,na.action=na.exclude)

where "dv" is the values from one raster and "iv1", "iv2" "iv3" ... are values from other rasters (up to 10 variables)  with the same extent and resolution.  It seems I should be able to do this out-of-memory using the Raster package, but I am confused how.  Whether I create a brick, stack, or set of individual Raster objects, I cannot figure out how to send the variables to the lm function without using getValues and thus calling everything into memory (mine cannot even handle two variables).

A point in the right direction would be much appreciated!


texts: How to calculate extreme point (most north or east etc) of geography(polygon) type in PostGIS

texts:
I have a table which include geography(polygon) type.

CREATE TABLE place (
  name varchar(100),
  geom GEOMETRY(POLYGON, 4326)
);
insert into place
SELECT 'place1', ST_GeomFromText('POLYGON((-71.1776585052917
 42.3902909739571,-71.1776820268866 42.3903701743239,
-71.1776063012595 42.3903825660754,-71.1775826583081
42.3903033653531,-71.1776585052917 42.3902909739571))',4326);

I want to know following points of geometry(polygon)

most north point
most south point
most east point
most west point

( This seems to be commonly called "extreme point")

For example, in the above image, the coordinates (x, y) of the four corners of the rectangle are the information I want to know.
(But my applications data is not true rectangle.that has more corner.)

I checked document, I couldn't find suitable function.

https://postgis.net/docs/reference.html

Please tell me how to get or calculate extreme points using SQL.


texts: Information of two attribute tables together in one

texts:
I need the information of a layer with track points (X and Y coordinates) and the information of a polygon layer (area) together in one layer. The problem is that the two attribute tables have no column or description in common, but I need the polygons and the trackpoints which create the outline of the polygon, in one Excel sheet.

Is there a chance to bring those layers together in one?


texts:
Hoping someone can help me out with this issue.

I've got a table of polygons which are grid squares which get overlaid on a map. These are used when querying data so that points get grouped into each grid square. All this works fine till a point is on the boundary between the grid squares. I've included an example of this below. If I use ST_Contains I get no results as the point isn't contained within any of the polygons. If I use ST_Intersects the point intersects with the boundaries of multiple grid squares which isn't suitable.

Is there a way I can ensure that the point is only considered to be in one polygon? In are previous system before migrating to using PostGIS we used the PNPOLY algorithm that ensured this.

DROP TABLE IF EXISTS areasTemp;
    CREATE TEMP TABLE areasTemp
    (
      id integer,
      areacode character varying(10),
      polygon geometry
    );
    INSERT INTO areasTemp VALUES (634,'1x50',ST_GeomFromText('POLYGON((1 50,2 50,2 51,1 51,1 50))', 4326)),
    (663,'0x51',ST_GeomFromText('POLYGON((0 51,1 51,1 52,0 52,0 51))', 4326)),
    (633,'1x51',ST_GeomFromText('POLYGON((1 51,2 51,2 52,1 52,1 51))', 4326)),
    (664,'0x50',ST_GeomFromText('POLYGON((0 50,1 50,1 51,0 51,0 50))', 4326));

-

SELECT * FROM areasTemp WHERE ST_Contains(ST_GeomFromText('POINT(1 51)', 4326), polygon);

0 rows retrieved

SELECT * FROM areasTemp WHERE ST_Intersects(ST_GeomFromText('POINT(1 51)', 4326), polygon);

4 rows retrieved


texts: Create multiple layers from one table on GeoServer?

texts:
I have a geometric SQL Server table with the following structure:

|--------|----------|---------|
|   id   |   geom   |  year   |
|--------|----------|---------|

I need to create 3 different layers in GeoServer from this table, filtering each by the year field, so that:

Layer A: contains the points with year 2000
Layer B: contains the points with year 2010
Layer C: contains the points with year 2020

This could be an example:

|--------|----------|---------|
|   id   |   geom   |  year   |
|--------|----------|---------|
|   10   |   foo    |   2000  |  ---> Layer A
|--------|----------|---------|
|   11   |   foo    |   2010  |  ---> Layer B
|--------|----------|---------|
|   12   |   foo    |   2010  |  ---> Layer B
|--------|----------|---------|
|   13   |   foo    |   2020  |  ---> Layer C
|--------|----------|---------|
|   14   |   foo    |   2020  |  ---> Layer C
|--------|----------|---------|

Is it possible to do this directly in GeoServer without creating 3
views in the database and importing them as independent layers?


texts: Problem creating TIN from multiple contour shapefiles with ArcGIS 9.3 Spatial Analyst

texts:
I have about 20 shapefiles that contain contours lines, and I'm trying to convert them to one TIN using Spatial Analyst in ArcGIS 9.3. I included all of the SHP's as input to generate a TIN, but it didn't work. I have tried with hard line, soft line, and mass points options, but nothing is working.

I think if I first created a single contour shapefile from the 20 then I would be able to generate a TIN. My question is: how can I create 1 contour layer from the original 20, or how can I create one TIN from multiple shapefiles?


texts: Selecting features from partial string in where clause with variable of Select Layer By Attribute in ArcPy

texts:
I am new to coding using ArcPy.
My code is:
roadname = "Ly Thuong Kiet"
arcpy.SelectLayerByAttribute_management("inLayer_lyr","NEW_SELECTION",""" "DIA_CHI" LIKE '{0}' """.format(roadname))

This code can run easily. However, I want to choose all the elements belonging to "DIA_CHI" field that have a name in the roadname.
Example: "Hem 32 Ly Thuong Kiet" will be contained in ("Ly Thuong Kiet").
When I wrote this code:
arcpy.SelectLayerByAttribute_management("inLayer_lyr","NEW_SELECTION",""" "DIA_CHI" LIKE "%Ly Thuong Kiet")

It can run. But when I use variable , I couldn't use percentage sign (%) in format.


texts: Point within, query in a buffer polygon in WMS like in ArcGIS

texts:
I need to do a query to filter the result inside a polygon or geometry, like this:
http://help.arcgis.com/en/webapi/javascript/arcgis/jssamples/query_buffer.html

But, I need to do this in a OpenGeo WMS server.

I have the following query:

  suite.opengeo.org/geoserver/usa/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&LAYERS=usa%3Astates&QUERY_LAYERS=usa%3Astates&STYLES=&BBOX=-79.89257812499999%2C41.68111756290652%2C-64.88525390624999%2C46.41513877649202&FEATURE_COUNT=5&HEIGHT=600&WIDTH=1366&FORMAT=image%2Fpng&INFO_FORMAT=text%2Fhtml&SRS=EPSG%3A4326&X=681&Y=298

And I need to filter the results inside of a polygon.
The parsed query:

  SERVICE: "WMS",
  VERSION: "1.1.1",
  REQUEST: "GetFeatureInfo",
  LAYERS: "usa:states",
  QUERY_LAYERS: "usa:states",
  STYLES: "",
  BBOX: BBOX,
  FEATURE_COUNT: "5",
  HEIGHT: HEIGHT,
  WIDTH: WIDTH,
  FORMAT: "image/png",
  INFO_FORMAT: "text/html",
  SRS: "EPSG:4326",
  X: X,
  Y: Y

I find this XML code, but I don't know if is this way:

FILTER: '<Filter xmlns:gml="http://www.opengis.net/gml"><Intersects><PropertyName>the_geom</PropertyName><gml:Polygon srsName="http://www.opengis.net/gml/srs/epsg.xml#2400"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>1 1</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></Intersects></Filter>'

EDIT:

In ESRI WMS, we can use this:

**point_within**
Description: An SLD style that uses a "within" spatial filter to render only features within a polygon geometry
URL of SLD file: http://sampleserver1.arcgisonline.com/arcgis/wms/slds/point_within.xml
GetMap request: http://sampleserver1.arcgisonline.com/ArcGIS/services/Specialty/ESRI_StatesCitiesRivers_USA/MapServer/WMSServer?VERSION=1.3.0&REQUEST=GetMap&CRS=CRS:84&BBOX=-178.217598,18.924782,-66.969271,71.406235&WIDTH=765&HEIGHT=360&LAYERS=0,1,2&STYLES=,,within&EXCEPTIONS=application/vnd.ogc.se_xml&FORMAT=image/png&BGCOLOR=0xFFFFFF&TRANSPARENT=TRUE&SLD=http://sampleserver1.arcgisonline.com/arcgis/wms/slds/point_within.xml

But, how I can do the same with others WMS?

By example the http://suite.opengeo.org/geoserver/usa/wms?

I tried using a FILTER with GML, but I dont know how to do this.
Someone help me?


texts:
I have Himawari-8 data (.nc files)
but I want to know the information such as projection, calibration, data acquisition,
I got this data from Meteorology agency in my country and I thought if the data has been modified so the original header is replaced with the new one (after the processing)
How do i know or read the original header?


texts: What exactly does hill shade or Shaded relief , if performed on SRTM DEM will output?

texts:
Two questions actually.
Does shaded relief or hill shade outputs mean "expelling of shadows"?

For a vast country like India, how to prepare a map from SRTM DEM, which does not have any shadow effects?


texts: Is there a programmatic way to swap attributes on a feature?

texts:
I am still pretty new to python (and programming for that matter) but I believe there must be an easier way to do a certain task. My task is to flip road centerline directions. The centerlines have address ranges in the attributes, and I need to flip them, from left to right. Basically I am swapping the value of some attributes.

Is that something I could accomplish with some scripting? I have a Python book on my desk that I've been reading through. I'm eager to tackle new problems with code.


texts: How to visuilze PostGIS data in GeoJSON format with PHP

texts:
I need to visualize PostGIS Geometry column as an overlay of a Leaflet tiled Map. I started a question here ,but didn't get any simple and satisfied answer. So I am starting this question. I got an option on Github to convert the PostGIS column to GeoJSON which can be visible through leaflet map.

I used the url as follows

http://localhost:8080/pgsql/PostGIS_GeoJSON.php?geotable=bgd_adm3&geomfield=geom

but it returns the following errors

Notice: Undefined index: parameters in C:\wamp\www\pgsql\PostGIS_GeoJSON.php on line 52

Notice: Undefined index: orderby in C:\wamp\www\pgsql\PostGIS_GeoJSON.php on line 54

Notice: Undefined index: limit in C:\wamp\www\pgsql\PostGIS_GeoJSON.php on line 61

Notice: Undefined index: offset in C:\wamp\www\pgsql\PostGIS_GeoJSON.php on line 63

Warning: pg_query() [function.pg-query]: Query failed: ERROR: function transform(geometry, integer) does not exist LINE 1: SELECT *, st_asgeojson(transform(geom,4326)) AS geojson FROM... ^ HINT: No function matches the given name and argument types. You might need to add explicit type casts. in C:\wamp\www\pgsql\PostGIS_GeoJSON.php on line 89
An SQL error occured.

What I have to do?


texts:
I'm making some test for the geocoder comparing it with google and I have noticed that the results depend a lot on the address format.

E.G.

These give addresses:

358 SEVENTH AVENUE MANHATTAN 10001
955 2ND AVENUE MANHATTAN 10022
530 AUDUBON AVENUE MANHATTAN 10040

Gives a wrong results, while these work just fine:

358 7th Ave, New York, NY 10001
955 2nd Ave, New York, NY 10022
530 Audubon Ave, New York, NY 10040

How reliable will be to go from the first group the second one using a library like this one and some strings mapping when required?


texts: Zonal statistics using Python console in QGIS gives AttributeError?

texts:
Here is how the log shows for one raster image. I want to do this same thing for several rasters. The batch processing is not working.

QGIS version: 3.8.0-Zanzibar
QGIS code revision: 11aff65f10
Qt version: 5.11.2
GDAL version: 2.4.1
GEOS version: 3.7.2-CAPI-1.11.0 b55d2125
PROJ version: Rel. 5.2.0, September 15th, 2018
Processing algorithm…
Algorithm 'Zonal statistics' starting…
Input parameters:
{ 'COLUMN_PREFIX' : 'sum_201204', 'INPUT_RASTER' : 'C:/Users/ghosh/Downloads/vbd_China_fishingban_analysis/tif/VBD_npp_20120401-20120430_asia_qf1-2-3-8-10-pc_v23_c201806110200.n_detect.china_crop.tif', 'INPUT_VECTOR' : 'C:/Users/ghosh/Downloads/vbd_China_fishingban_analysis/vectors/china_eez_v9_banarea.shp', 'RASTER_BAND' : 1, 'STATS' : [1] }

Execution completed in 17.04 seconds
Results:
{'INPUT_VECTOR': <qgis._core.QgsVectorLayer object at 0x0000029580A2E4C8>}

Loading resulting layers
Algorithm 'Zonal statistics' finished

This is the error for batch processing

  Input parameters: {} Traceback (most recent call last):

  File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\qgis\ZonalStatistics.py", line 128, in prepareAlgorithm self.raster_interface = rasterLayer.dataProvider().clone() AttributeError: 'NoneType' object has no attribute 'dataProvider'

  Batch execution completed in 0.02 seconds


texts: Converting Excel files into geodatabase with ArcPy?

texts:
I have several Excel files that contain several sheets. I need to convert these files into one table in gdb. I wrote script that convert the whole sheet in separate table. How can I export a range (f.e. A8:G(rows on the sheet)) of the sheet to one table.

import arcpy, os, xlrd,xlwt
from arcpy import env

arcpy.env.overwriteOutput = True

fname = 'D:/excel_files/exc4.xls'
output_gdb = "D:/P01_ExctoGDB/GDB_EXC.gdb"
# Open the workbook
rb = xlrd.open_workbook(fname)

# List sheet names, and pull a sheet by name
sheet_names = rb.sheet_names()

for i in range(0,len(sheet_names)):
    sheet=rb.sheet_by_index(i)
    val=sheet.row_values(0)[0]
    arcpy.AddMessage(val)
    vals=[sheet.row_values(8) for rownum in range(sheet.nrows)]

    outTablename= output_gdb+"/"+str(sheet_names[i]).replace(" ",'').replace(".","")
    arcpy.ExcelToTable_conversion(fname,outTablename, sheet_names[i])


texts: SQL Expression with Multiple Values Parameter in ArcPy

texts:
I`m trying to select many values from the parameter "Estado" which will be use in the query in the code. This parameter has a value list, string type.

I think it`s a format() problem...But also I will have problems to put in the SQL expression format: 'VALUE1','VALUE2'.

I tried with .split(';'), but do not work.

import arcpy

Layer1 = arcpy.GetParameterAsText(0)
Tabela1 = arcpy.GetParameterAsText(1)
Estado = arcpy.GetParameterAsText(2).split(';')
Ano = arcpy.GetParameterAsText(3).split(';')
mxd = arcpy.mapping.MapDocument("CURRENT")
Saida = arcpy.GetParameterAsText(4)
stats = []

for lyr in arcpy.mapping.ListLayers(mxd):
    lyr.name = Layer1
    lyr.definitionQuery = """UF IN ({})""".format(Estado)

for field in arcpy.ListFields(Tabela1):
    if field.name in (Ano):
        stats.append([field.name, "SUM"])

arcpy.Statistics_analysis(Layer1, Saida, stats, "")
arcpy.RefreshActiveView()
arcpy.RefreshTOC()

When I tried to run the script, I got this error:

Traceback (most recent call last):
  File "F:\AgroBD.gdb\M6a.py", line 19, in <module>
    arcpy.Statistics_analysis(Layer1, Saida, stats, "")
 File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\analysis.py",
line 1160, in Statistics
    raise e
ExecuteError: ERROR 999999: Error executing function.
An invalid SQL statement was used.
An invalid SQL statement was used. [Brasil_Municipios]
An invalid SQL statement was used. [SELECT * FROM Brasil_Municipios WHERE UF
IN ([u'AC', u'AM'])]
An invalid SQL statement was used.
An invalid SQL statement was used. [Brasil_Municipios]
An invalid SQL statement was used. [SELECT OBJECTID FROM Brasil_Municipios
WHERE UF IN ([u'AC', u'AM'])]
The table was not found. [out3]
Failed to execute (Statistics).

Failed to execute (Script1).


texts: How to dissolve overlapping distance isochrones

texts:
For my master's thesis I need to create distance isochrones (10/20/30/../100km) around certain coordinates. In my case they overlap and I want the overlapping parts to dissolve, but at the same time keep the distance intervals of 10km. Is there a smooth way to do this?
Here is the example pic of the overlapping isochrones:

You can see there is an underlying isochrone with the according distances.. is there someway to merge them and still get the according distance intervals?


texts: Add Action (-Button) to Attributes Form (QGIS)

texts:
Is it possible, to simply add a Action to a Attribute Form (Drag'n Drop Designer)? I know, I could use the QT-Designer and Python, but I'm looking for a fast and simple way without python knowledge.

And to get further: the Action should be used to delete the current feature.


texts: Specifying Resource Field For Labeling in ArcPy with ArcMap

texts:
How do I specify which field acts as resource for labeling a layer using ArcPy?

Following code is some settings that I am doing on a standalone ArcPy snippet and as you can see I added layer.showLabels = True but assuming I have a field named name I do not know how to add it to the code.

layer = arcpy.mapping.Layer("Primary")
layer.minScale = 10000
layer.showLabels = True


texts:
I have a client, that needs some geological maps digitized (polygons) in a Cartesian grid.  It would be a basic, square grid with North up, and coordinates in x and y.  Some of the maps will need to be digitized in meters and some in feet.  I've looked around a lot for a way to set this up in QGIS, but I can't find anything that I can understand.

I think if my technical knowledge was better, I could make a custom projection that would work.  I'm baffled that this answer isn't easier to find.  Although I wouldn't be surprised if I'm overlooking something.

I'm not interested in placing this spatially.  It can stay in the Cartesian grid.  I'm doing the digitizing in QGIS and he is placing the SHP files into MapInfo for plotting.

Edit: A question that I've been wondering is if it is fine to digitize in a UTM grid as a SHP file, then just delete the PRJ?  Would the SHP file then have the right coordinates in a flat plane?  Or is there more to it than that?

Thanks,


texts: Field Calculator: extract single cell based on it's OBJECTID?

texts:
Is it possible to copy a single value inside a field based on it's OBJECTID from another field?

Some sort of if-then-else logic?

For example I would like to copy the value in the Area field corresponding to OBJECTID 3. And to have it in the first row in the Level field with the rest of it's rows remaining to 'Null'.

I am using ArcGIS Pro 1.4 with Python 3.5.


texts: How to draw line/multiline, polygon/multipolygon with SharpMap from Coordinate Source

texts:
I have a list of n points, exactly x objects in GeoJSON from a CouchDB. I want to generate lines/multilines and polygons/mutipolygons with
SharpMap and add this to a Layer in VB.NET. For points its easy, but for the rest: I don't get it. Must I use ICoordinateSequence?

So, can anybody post here some example to generate such geometries dynamically from coord-list-like text- or GeoJSON/JSON-Sources (maybe here its another workaround than for plain text?), much detail as possible.


texts: Acquiring Digital Elevation Models for Argentina from commercial providers

texts:
I'm interested in obtaining Digital Elevation Models for areas in the 20,000 acre range, as I use them to create contour maps for agriculture.

The USGS Earth Explorer has a wonderful free DEM database with world coverage, which is important to me, since I live in Argentina. Sadly, the resolution for the global dataset is too coarse, so it's useful to analyze whole countries, but not for smaller scale analysis.

So I searched the web looking for commercial satellite data providers that offer DEM files.

I found two companies that offer satellite imagery and derived products, like preprocessed DEMs, with worldwide coverage:

Airbus Defense & Space (http://www.geo-airbusds.com/worlddem-sampledata/)
BlackBridge ( http://blackbridge.com/rapideye/all-products/nextmap.html )

Of course, the data is not cheap, and, at least in the case of Airbus, they'll not take orders of less than 500 km2.

I wonder if there are any other providers of this type of data, ideally free and open, but realistically, commercial, that may be better tailored for smaller scale projects. A resolution of around 12m is perfect, but in the case of free/open sources, of course beggars can't be choosers :)


texts: Querying featureLayer using date field and ArcGIS API for JavaScript?

texts:
I can query a featureLayer in my javascript map to only show features meeting criteria for normal attributes (strings/integers), as below for example

featurelayer.setDefinitionExpression("status = 0"); // working fine

However i would like to query based on date field in my feature layer -  (e.g. only display records newer than 01/01/2014)

in the service metadata the field is described as:

date_created ( type: esriFieldTypeDate , alias: DATE_CREATED , length: 36 )

I have looked through the the documentation, however i have seen no mention of how to do this kind of query, the closest thing mentioned is the setTimeDefinition, but this appears more related to making time-aware maps rather than performing a simple filter on a layer


texts: Connect a Google Spreadsheet for importing a Feature Collection to Google Earth Engine

texts: I have a google spreadsheet consisting a feature collection table. Is it possible to import the table in Google Earth Engine? I am a new learner of Earth Engine and was thinking that we have to upload a shapefile to the "Assets" and then import it to the code editor. So, can we just import it directly from a google spreadsheet using the url?


texts:
I'm new to PostGIS (and SQL in general), and am running through the foss4g workshop, using their nyc data.
My code below works as intended : For every metro stop, I calculate the number of streets in a 100m radius and 200m radius.

CREATE TABLE close200 AS SELECT
nyc_subway_stations.name,
     count(DISTINCT nyc_streets.name) count200
FROM nyc_subway_stations, nyc_streets
WHERE  ST_DWITHIN(nyc_subway_stations.geom, nyc_streets.geom, 200)
GROUP BY nyc_subway_stations.name
ORDER BY count200 DESC;

CREATE TABLE close100 AS SELECT
nyc_subway_stations.name,
     count(DISTINCT nyc_streets.name) count100
FROM nyc_subway_stations, nyc_streets
WHERE  ST_DWITHIN(nyc_subway_stations.geom, nyc_streets.geom, 100)
GROUP BY nyc_subway_stations.name
ORDER BY count100 DESC;

CREATE TABLE proxMetro AS SELECT
close100.name, close100.count100, close200.count200
FROM close100, close200
WHERE close100.name = close200.name;
DROP TABLE close100, close200;
SELECT * FROM proxMetro;

I feel that this isn't the best way to write this out, but I'm not sure how to shorten my code or make it less heavy. I shouldn't have to create two extra tables and should be able to select the information directly from my original tables. My very first intuition was to use nested queries :

SELECT
    nyc_subway_stations.name,
    (SELECT count(DISTINCT nyc_streets.name) count100
        FROM nyc_subway_stations, nyc_streets
        WHERE  ST_DWITHIN(nyc_subway_stations.geom, nyc_streets.geom, 100)),
    (SELECT count(DISTINCT nyc_streets.name) count200
        FROM nyc_subway_stations, nyc_streets
        WHERE  ST_DWITHIN(nyc_subway_stations.geom, nyc_streets.geom, 200))
FROM nyc_subway_stations, nyc_streets
GROUP BY nyc_subway_stations.name
ORDER BY count100 DESC;

Which returns a same number for every row, which I'm not sure corresponds to anything. Is this second method on the right track? If so what am I doing wrong?

I'm using PostGreSQL 9 with PostGis 2


texts: Setting the name of feature in geoJSON and retrieving client side

texts:
How do I extract the name of a feature? I want to make it so when I mouse over the feature it is highlighted and is able to have a tooltip with the name of the feature that is moused over.

Here is part of my geoJSON;

"features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -13623642.663851619,
              5984966.7642723517
            ],
            [
              -13623264.36027813,
              5984616.0795468809
            ],
            [
              -13623234.20348176,
              5984616.8260022318
            ],
            [
              -13623255.850687083,
              5984636.9802967971
            ]
          ]
        ]
      },
      "name": "testName"
    },
    },
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [

And the code I am sing to try to get it;

var selectPointerMove = new ol.interaction.Select({
                            condition: ol.events.condition.pointerMove
                        });
                        map.addInteraction(selectPointerMove);
                        selectPointerMove.on('select', function (e) {
                            console.log(e.target.getFeatures().getLength())
                            console.log(e.target.getFeatures().get('name'))
                        });

Seems like whatever I try I cannot extract the name.

I resolved the issue by adding the properties to my geoJSON so it looks like this;

    "features": [
    {
      "type": "Feature",
      "properties": {
        "name": "IMPORT"
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -13623642.663851619,
              5984966.7642723517
            ],
            [
              -13623264.36027813,
              5984616.0795468809
            ],
            [
              -13623234.20348176,
              5984616.8260022318
            ],
            [
              -13623255.850687083,
              5984636.9802967971
            ]
          ]
        ]
      }
    },
    },
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [

And using the Javascript;

map.on('pointermove', function (evt) {
    if (evt.dragging) {
        return;
    }
    pixel = map.getEventPixel(evt.originalEvent);
});
var feature = map.forEachFeatureAtPixel(pixel, function (feature) {
        return feature;
    });
console.log(feature.get('name'))


texts: Merge without dissolve in ArcGIS for Desktop?

texts: I am trying to merge several polygons in a feature dataset into one. Normally I would use a simple merge but the issue is that merge automatically dissolves. I on the other hand want to have them all represented by one single record, but maintain their existing boundaries.The top image is what happens. But the bottom image is what I want. As mentioned originally the left would represent by two separate selectable polygons and the right is the output being one record/one selectable polygon.


texts: How to find within a geodatabase from the catalog window

texts:
I have a geodatabase containing hundreds of feature classes that I want to be able to search within so that I know which feature class to load.  All the feature classes come from different data sources, so there is no consistency in naming or structure.

As an example, from my hundred or so feature classes within my geodatabase, say I have two feature classes structured like this:

FC1

ID Animal Age Sex

01 Cat    3   F

02 Dog    3   M

FC2

ID Legs Genus Age

01 3    Dog   16

02 6.5  Arach 3

What I want is to search for Dog within the catalog window, and get a list of feature classes that have the word Dog in them, and from there I can work out which feature classes to load.

Is this possible in Arc?


texts:
A collaborator has given me coordinates of his sites in NAD27 Mexico. How do I go about converting it to Lat Long? I don't have ArcGIS, but is there a way to do the conversion in that program?

Edit: I have looked into a couple places that did State Plane coordinates, but did not find ones specifically for Mexico. I will try qgis.


texts: Handling dangling nodes in ArcGIS Desktop?

texts:
I have a datafile of line features. The creation of a topology revealed that there are too many dangling nodes to be corrected manually.

I was thinking of snapping to the closest line feature in an automated process, however, I do not have the Data Reviewer license of ArcGIS with which this would actually be possible (with the tool 'Snap').

How do I correct those dangling nodes in an automated process probably by a work around without the Data Reviewer license?


texts: Determining driver, spatial reference and other props from an open GDAL layer

texts:
I have a very large network share full of data (raster, vector, sometimes mixed inside of geopackage and/or file geodatabase containers).  I have a set of classes that crawl through it all and use arcpy.Describe() to extract from each dataset its datatype, extent and spatial reference system.  That then goes into a single enterprise db layer that shows our spatial holdings and a bit of heritage data of each element.
I'd like to try making it open source for other users to run who do not have licenses. GDAL also offers me some more insight into what specific driver is being used, to further categorize the data, theoretically.
I'd like to avoid hardcoding around file extensions and ask GDAL (or fiona, or shapely, not picky here) instead - [update] because I don't know what sort of files the indexer will encounter.  Assuming a subset of the data, I can do something like this:
import gdal

filelist = [
    '/share/Data/Terrain/N43W072.hgt',
    '/share/Data/Terrain/N43W072_StrahlerOrder_8.tif',
    '/share/Data/Terrain/N43W072_Channels_Clipped.shp',
    '/share/Data/Terrain/N43W072_Channels.shp',
    '/share/Data/Network/NH/nh_road_network.gdb',
    '/share/Data/Network/MA/ma_road_network.gpkg'
]

for f in filelist:
    try:
        dataset = gdal.Open(f)
    except:
        # do stuff

but I am not clear how to then summon the properties of dataset (I'm presuming that the open dataset object holds the elements I'm looking for). I've looked at https://gdal.org/python/ but there is nothing I can find that says "here is how you get all the properties of an gdal.Open() object."
Ultimately the properties I want to see would be presented something like this:
{
    "driver": "GPKG",
    "driver_long_name": "GeoPackage raster"
    "type": "Raster",
    "extent": [(miny, minx), (maxy, minx), (maxy, maxx), (miny, maxx)]
    "CRS" : {'epsg:4326'}
}


texts: I have some GPX data for flights I've taken, but the waypoints are very intermittent - often times 100's of km apart.  I have an app called Fog of the World which I use as a hobby to track the places I've been.  It allows me to import GPX data, but will not create a line when there is more than 5km between waypoints.  Most GPX visualizers will create this line regardless.  I was wondering if there was an easy way to add extra waypoints between the existing way points (that follow the lines) to satisfy the needs of the app.  I've tried the hand drawing tools on some websites, but obviously that's such a tedious way to map out these flights.


texts: Join a table to a parcel layer and get duplicate parcels, reflecting individual leases

texts:
I have a table of 3,000+ lease files that I'm trying to join to our parcel layer. Each parcel has 4+ leases associated with it, and I need each to be represented by an actual polygon.  When I do a traditional join it joins 1 to 1, so it shows 1 of the 4 leases and the others do not join at all. Is it possible to join 1 parcel to many leases and have each show up as an individual polygon?
I'm using Arc 10.1.

Thanks in advance,
Matt


texts: Eclipse PyDev has quit recognizing QGIS objects

texts:
I followed the instructions for setting up Eclipse with PyDev to start creating QGIS plug ins found here: http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/ide_debugging.html

after doing that I was able to create a QGIS plugin and it works as intended, but when I opened the project this morning eclipse no longer recognizes the QgsPoint or QgsRaster objects, and is now putting it's standard red line under them to I have checked

    Window > Preferences > PyDev> Interpreters > Python

and made sure that Libraries still contains the QGIS and QGIS plugin references and that Force Builtins still contains qgis

This problem does not seem to be altering my ability to modify the code. I can change the code, and after a reload QGIS will run the plugin with the changes I make.

I would like to know how to get eclipse+PyDev to recognize the QGIS objects again so that I can get rid of the red lines and focus on actual errors while debugging.

Update:

The plot thickens. This screen shot was taken with my mouse hovering over the pixelData0 object and it clearly shows that Eclipse can resolve it to a layers.dataProvider.identify object from  QGISRasterLayerCreation so Eclipse clearly can see some of the QGIS objects, but if you look at the same line the QgsPoint reference is being underlined in red. When I hover over it the error message says Undefined Variable: QgsPoint


texts: Reverse XY coordinates from csv to qgis?

texts:
I find it very difficult to import GPS data from CSV to QGIS and finally I have discovered that the first column in CSV table was Y coordinate and the 2nd was X coordinate.

In the picture you can see the settings for the file import

To check the real position in google maps I have to insert (first value, second value) so I believe that was the same in QGIS

Why is this happening? Is it regular?


texts: Using DrawFeature to draw a line between features

texts:
I have a DrawFeature control that I want use to draw a line between two features (on a separate layer from the lines). I use a Snapping control to snap the cursor to the features.

Any ideas on how can the drawing be limited so that it only allows to end the line on a feature?

Illustration (O represents a feature):

O---O allowed
O---  not allowed
 ---  not allowed


texts: Where does the copied data get placed when publishing a service to ArcGIS Server

texts:
I use ArcGIS for Server 10.1. What is not clear to me is, where does server save data when I upload it to it?
I mean the case, when data in the service is not referenced, but it is inside the service definition file.

Is it saved in some folder?
Or is it saved in some geodatabase included in the server? Where?


texts: Converting GeoTIFF to array of normalized (0 - 1 range) values using C# GDAL Nuget Package

texts:
I currently have a Unity3D tool that allows users to import different heightmap files into the Unity3d Game Engine, in order to produce Unity Terrain objects. The tool currently has limited support for GeoTIFF files, as I am manually reading them and that only works for some projection types. I would like to expand support to (hopefully) all GeoTIFF files by using the GDAL Nuget package from Tamas Szekeres.

I basically need to take a GeoTIFF file, and produce a heightmap from it, which is a float[,] array of normalized values between 0 and 1. Unity uses this information to create the Terrain object.

What is the best way to do this? I first attempted to use the DataSet.ReadRaster method, but that is crashing Unity with an error saying I am trying to read an invalid memory address. Even if I could get that method working, I am not sure if it is the correct route to take.

Here is what that looked like:

    double[] rasterData = new double[outputWidth * outputLength];
    ds.ReadRaster(0, 0, 65, 65, rasterData, 65, 65, ds.RasterCount, null, 0, 0, 0);

Most suggestions I see say to use the individual Raster Bands, which I did test and it does not cause an error (I did not actually try to use this data), however how do I handle different GeoTIFF files (using different projections and data formats) with this method? For example, I have many files with a single Raster Band and which store individual height data as floats, which seems easy enough to read, but then I have other files where the height data is in RGB format.

Is there a single method I can call to translate the data, no matter what format it is, to a normalized value or a value that I can manually normalize? Or do I need to handle each projection type independently?

I am using this in a 64 bit Unity Editor.

I have been using LibTiff.Net to manually read the GeoTIFF files. I thought using GDAL I would be able to gain access to some sort of auto translate function to convert almost any GeoTIFF file to either an array of raw height values or even an array of normalized values. But the more I look into it it seems I will have to read the DataType of each file and manually convert each data point as I am already doing, so I am not sure if there is any advantage in switching to GDAL.

Is that correct?


texts: Convert from Alt/Az coordinates w.r.t. point on surface to (x,y,z) w.r.t Earth centre?

texts:
I know how an antenna located on Earth surface is pointing to an object in space, by means of Altitude and Azimuth.
I would like to show the pointed object in the 3d view, so, as far as I can understand, I need its position in (x,y,z) format (Cartesian3 object?)
But I can't understand exactly which is the exact process to follow.
I am using this example using CesiumJS library to make experiments; I added one point on surface by adding these lines inside function addMultiplePoints() :
viewer.entities.add({
position: Cesium.Cartesian3.fromDegrees(290.6, -35.78),
point: {
  color: Cesium.Color.RED,
  pixelSize: 18,
},
label: {
  text: "MLG",
  font: "14px Helvetica",
  pixelOffset: new Cesium.Cartesian2(0.0, 20),
},
description : "MLG"
});

Now I have to add a point with specific location w.r.t this point, say Alt=30, Az=0.
I am trying with this... but nothing happens; indeed, in this code there is no origin point for the coordinates...
  H2pointing = {"alt" : 30, "az" : 0};
  clock = 90 - H2pointing.az;
  cone = 90 - H2pointing.alt;
  H2spherical = new Cesium.Spherical(clock, cone, 10);
  console.log(clock);
  console.log(cone);
  console.log(H2spherical);
  H2cartesian = Cesium.Cartesian3.fromSpherical(H2spherical, null);

viewer.entities.add({
position: H2cartesian,
point: {
  color: Cesium.Color.BLUE,
  pixelSize: 18,
},
label: {
  text: "Hayabusa2",
  font: "14px Helvetica",
  pixelOffset: new Cesium.Cartesian2(0.0, 20),
},
description : "Hayabusa2"
});

So, assuming my calculations are right (are they?), how do I center this point on previous point?
Maybe I could use another method:
Cesium.Transforms.headingPitchRollToFixedFrame(origin, headingPitchRoll, ellipsoid, fixedFrameTransform, result)

If I specify Roll = 0, I should get the equivalent of Altitude and Azimuth.
But, also in this case, how to apply this transform to a point?

headingPitchRollToFixedFrame Computes a 4x4 transformation matrix from a reference frame with axes
computed from the heading-pitch-roll angles centered at the provided
origin to the provided ellipsoid's fixed reference frame. Heading is
the rotation from the local north direction where a positive angle is
increasing eastward. Pitch is the rotation from the local east-north
plane. Positive pitch angles are above the plane. Negative pitch
angles are below the plane. Roll is the first rotation applied about
the local east axis.


texts: How to calculate max drift along converging lines?

texts:
I have a plot of 2 navigations, red and blue, of two ships that are converging and diverging.  Most of the time they are almost parallel.
Is there a way to calculate the maximum drift along the 2 lines?

NOTE -  the command must take into account the drift over a certain range  only,  say on a 100m range, or else, as there are different passages of the lines, enter image description here the command would see also the farmost drift from lines too far apart... please see the photos