texts:
I want to overlay some data whose projection is WGS-84 on Google map layer in OpenLayers. But I just can't make them in the right place. I did as follows:

map = new OpenLayers.Map('map', {
        numZoomLevels: 20,
        projection: new OpenLayers.Projection("EPSG:900913"),
        displayProjection: new OpenLayers.Projection("EPSG: 4326")
        });
googlelayer = new OpenLayers.Layer.Google("Google street", {sphericalMercator: true});
map.addLayer(googlelayer);
veclayer = new OpenLayers.Layer.Vector("vector", {
                                    projection: map.displayProjection
                                    };
var geojson_format = new OpenLayers.Format.GeoJSON();
veclayer.addFeatures(geojson_format.read(jsonData));

Though I have assigned veclayer in 4326 projection, but it is still interpreted as 900913, and the display coordination system is also 900913, though I set displayProjection to 4326. What mistake do I make?


texts: ArcGIS Service Backend Deletion: does anyone have instructions to remove a service from the server itself (filesystem/registry/etc)?

texts:
I have a rogue service on a a single-server ArcGIS Enterprise instance which, after spending quite a bit of time on the horn with ESRI, I've been unable to delete (or stop, for that matter).  We've looked at locks, missing files, and permissions, and of course tried everything under the sun from the Server Manager, Server Administrator, and ArcGIS Catalog interfaces.

Surely, there must be a workflow to clear out an angry service from the machine itself.


texts: Cutting off polygons outside of border in Leaflet?

texts:
I have an application where I create polygons based on values I calculate from points on a map. The polygons are then overlayed on top of a city using a single leaflet FeatureGroup. Another one is then used to draw the border according to a GeoJSON I have with borders for each of the cities where I create these polygons. Also, each polygon belongs to a tier in a tiered system, much like an isopleth.

However, I'd like to cut these polygons so that only what is located inside the city borders is shown. Are there any tricks using Leaflet to achieve that? If at all possible, I'd like to avoid having to use any other package and need to adapt the way my polygons are defined to work with this other package.

This is how my application is looking at the moment:

Since I was using a python server to generate my polygons, I ended up using Shapely to calculate the polygon intersections and only show what's inside the borders.


texts: Custom Attribute Display Pop-up (hyperlink converting "/' to "%2F")

texts:
I have created a custom attribute display links in an ArcGIS Online webmap that include a field from my ArcGIS Server Feature Service attribute table (using attribute-based URL parameters). I am having an issue where the "/" in my attribute field record is appearing as "%2F" in the hyperlink and thus the proper URL is not opening.
Any suggestions on this?


texts: Is it possible to create vector style in QGIS that is a border inside border?

texts:
In England, when submitting applications for Countryside Stewardship Scheme woodland grants you are required to mark on a paper map, supplied by the authority, the woodlands that are to be included.

The requirements are very specific, they require a red outside boundary, with a yellow inside, as in the guidance notes that come with the maps. See the image below, where I have added the colour with felt tip pens.

For large estates with many scattered woodlands this can be very time consuming, and if you make a mistake in your colouring in, you really have to start all over again.

As I already have all the areas as polygons in vector shapefiles I have scanned the maps and georeferenced them into my QGIS project.

Is it possible to replicate this border style of red on the outside and yellow inside that?

I have tried adding another line style for the yellow, and making it wider, but of course this extends both sides of the red border.

I have found a similar question for ArcGIS Using directional shading on boundaries for cartographic effect in ArcGIS Desktop?


texts: How to efficiently compute if point is on polyline?

texts:
I have a CSV file with data in GeoJSON format with 10k+ lines with various open / public data I want to use. So far it's in a file but later it may be in a SQLite database or even better be kept in RAM. This is to be used in a mobile app.

There are many columns in that file but what interests me are the geo shape column which is made of polylines (linestring with several points) and the column called "geo center"(it is has the center between start and end points of the polyline).

I want to know if a location (provided by the mobile device) is on any polylines defined in the file.

So my first thought is to filter all geo centers that are near enough to the given location, and then for each filtered polylines, check if the given location belongs the series of point couples.

Eg : location is L

Polyline 1 has 3 geo points X Y Z

My algorithm is (after filtering step) :

Is L on [X Y] ?
If not is L on [Y Z] ?
...

Is it the most efficient way to tell if a point is on a polyline or is there any other way ?


texts: Making density map based on number of hospitals in statistical area using ArcGIS Pro

texts:
I have the longitude and latitude of hospitals in Australia and I want to make a density plot (Symbology - graduated colours) to visualise the distribution of hospitals by a statistical area (defined by Australian Bureau of Stats).
I've downloaded the shapefile from Aus Bureau of Stats.
The hospital data looks like this and I import it as XY event data:
Hospital name,Street address,Suburb,Postcode,State,Local Hospital Network (LHN),Primary Health Network area (PHN),Website,Description,Sector,Beds,Latitude,Longitude
Abbotsford Private Hospital,61 Cambridge Street,West Leederville,6007,WA,,,www.abbotsfordhospital.com.au,,Private,<50,-31.940992,115.837344
Adelaide Clinic,33 Park Terrace,Gilberton,5081,SA,,,www.adelaideclinic.com.au/,,Private,50-99,-34.90658,138.613167

I can get the hospitals overlayed on the map but I cant work out how to get it to colour the statistical area based on the number of hospitals in this area.


texts:
I am trying to convert a line shapefile with branches into just one polyline, actually deleting every other diramation.

Imagine something like this:

:________|__________|____|||______

into

:___________________________________

(dont consider the two dots, had to put them for format reason)

I could do it manually but it would take alot of time. Maybe using some sort pf shortest path could work but I am not sure.


texts: Using Qt Designer to Create a Form in QGIS that can search for a specific point in a layer

texts:
I would like to ask if you can create a form in Qt Designer to create a form that can search for a point on the map using data that the user will provide such as the Name.

I am currently creating a map in QGIS for my school project and i was thinking if i can create a form that can search the layer for a specific point and show the result in a dialog box using Qt Designer.


texts:
I was a little surprised to see the output of a .jpg file when gdalinfo said it had overviews.

Can jpegs really have overviews or what does this mean? Can they be removed to reduce file size?

Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
  SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 4912.0)
Upper Right ( 7360.0,    0.0)
Lower Right ( 7360.0, 4912.0)
Center      ( 3680.0, 2456.0)
Band 1 Block=7360x1 Type=Byte, ColorInterp=Red
  Overviews: 3680x2456, 1840x1228, 920x614
  Image Structure Metadata:
    COMPRESSION=JPEG
Band 2 Block=7360x1 Type=Byte, ColorInterp=Green
  Overviews: 3680x2456, 1840x1228, 920x614
  Image Structure Metadata:
    COMPRESSION=JPEG
Band 3 Block=7360x1 Type=Byte, ColorInterp=Blue
  Overviews: 3680x2456, 1840x1228, 920x614
  Image Structure Metadata:
    COMPRESSION=JPEG


texts: How to delete custom color ramp in QGIS?

texts: I have customized a new color ramp in qgis. Now I want to delete that. I could not find any way to delete this.


texts:
I am trying to install GDAL on the remote server.
I do:

pip install GDAL==1.9.1

It gets installed and then I try to do:

gdalinfo --version

But I get 'command not found'.

I understand that I have to write the environment variable in the .bash_profile script.

So I open the file with vim and I type:

export PATH=/Library/Frameworks/GDAL.framework/Programs:$PATH

And then:

source ~/.bash_profile

But still doesn't work.

I checked the location of my GDAL with: which ogr2ogr command and I get:

/usr/bin/which: **no ogr2ogr in** (/Library/Frameworks/GDAL.framework/Programs:/Library/Frameworks/Python.framework/Versions/2.7/bin:/Library/Frameworks/GDAL.framework/Programs:/Library/Frameworks/GDAL.framework/Programs:/Library/Frameworks/GDAL.framework/Programs:/home/vagrant/.venvs/geonode/bin:/usr/local/bin:/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/sbin:/home/vagrant/bin:/home/vagrant/bin:/home/vagrant/bin:/home/vagrant/bin:/home/vagrant/bin:/home/vagrant/bin:/home/vagrant/bin:/home/vagrant/bin)


texts:
I have a location off of the coast of Madagascar with longitude and latitude:
X = 43.28191425

Y = -22.49701687

This point also has a Z value (height). I want to convert the lat,lon to UTM co-ords and convert the Z value from ellipsoid height to Orthometric height (EGM2008)
I can convert the Z value using:
import pyproj
from pyproj import Proj
from pyproj import Transformer

transformerh = Transformer.from_crs("epsg:4326", "epsg:3855")
Xegm, Yegm, Zegm = transformerh.transform(X, Y, Z)

However, when I try converting the X and Y co-ords to UTM with:
    myProj = Proj("+proj=utm +zone=38, +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

    Xutm, Yutm = myProj(X, Y)

it gives me the co-ords:
Xutm = 323263.48408316
Yutm = -2488855.56162304

so the Yutm value is negative and makes no sense for its location, although the Xutm value is correct. The lat and lon point is definitely in UTM zone 38 (South) and when I try the conversion with other software (rsgsilib) it gives me the co-ords:
Xutm_check = 323263.4843173389
Yutm_check = 7511144.4386357535

which match the co-ords in Google Earth.
Am I incorrectly assigning myProj?


texts: Reading source code of QGIS plugins in Processsing Toolbox?

texts: I want to use the code snippets of "Generate points(pixel centroids) inside polygon" plugin in Processing Toolbox of QGIS. I tried searching in QGIS github, but couldn't find it.


texts: Modify spatial lines according to the minimum distance between line coordinates in R

texts:
I have an output of  getContourLines() from  contoureR package. The data frame output was modified and some points in the contour path were removed according to some masks, due to that: the distance between the points in the contour paths became longer as you can see in the leaflet map.
I used the function here to connect points to lines.
head of the data frame :
> head(output_df)
     LID GID PID Group        x          y z
249    1   1 220   1-1 165.0000  -4.010842 0
250    1   1 221   1-1 164.6938  -4.500000 0
251    1   1 222   1-1 164.3503  -5.350302 0
252    1   1 223   1-1 163.8953  -6.000000 0
1543   1  18 203  1-18 201.0000 -10.075164 0
1552   1  18 212  1-18 207.0000 -10.937481 0

I would like to separate any line that has a big gap between its points to smaller lines for example line 1-18 in the attached map, Constructing lines from x,y coordinations without ID or sorting didn't work.
Another approach is running loops over all points and filter lines to keep only lines coordinates that pass the threshold, this is done before converting to spatial lines. It removes too many points and can be run endlessly to reach the goal.
Loop code is :
#Filter points that are really far from their line group
end_of_loop<-nrow(output_df)-1
filter_of_far_points<-vector(length = end_of_loop)
for( j in 1:end_of_loop){

  distance_to_next_point <-sqrt((output_df$x[j]-output_df$x[j+1])^2+(output_df$y[j+1]-output_df$y[j+1])^2)
  if(output_df$Group[j]==output_df$Group[j+1]){
    if(distance_to_next_point< 4.5){
      filter_of_far_points[j]<-1
      #print(c(j,"yes"))
    }else{
      #print(c(j,"NO",distance_to_next_point))
      filter_of_far_points[j]<-0
    }
  }
}
#remove points which did not pass filter
nop<-which(filter_of_far_points==0)
output_df<-output_df[-nop,]

EDIT:
Can I transform it into Spatiallines and then cut lines to segments, where the distance between line coordinations exceeds the minimum distance threshold?


texts: Install PostGIS 2.0.3 on PostgreSQL x64 9.1 on Windows Server 2008 R2

texts: I just try to install the PostGIS 2.0.3 on PostgreSQL x64 9.1 on the Windows Server 2008 R2. I got the error message "Error openning file for writing" C:\Program Files\PostgreSQL\9.1\bin\libiconv-2.dll. Should I click "Abort", "Retry" or "Ignore"? Thanks. :)


texts: What algorithm allows me to calculate the set of all points inside a polygon?

texts:
Imagine I have a polygon, which represents the countours of a building on an integer 3D map (everything is a 1 meter by 1 meter block in that map). In the image below, you see several buildings without roofs.

In order to draw the roofs, I need to know what points are inside the polygon.

What algorithm can I use for that?

People mentioned convex hull, but it gives me the polygon, which I already have.

Update 1 (01.10.2015 13:34 MSK): I found following potential names:

Spatial joining (suggested by Joseph)
Scan line algorithm (another explanation)
Z-buffering

Update 2 (01.10.2015 14:22 MSK): Imagine I have a rectangle defined by the black line in the image below, drawn on a grid.

I need an algorithm, to find out the coordinates of all points inside the polygon, i. e. all the green points.

Update 3 (01.10.2015 15:09 MSK): Note that this is a Minecraft map, which you can imagine as a 3-dimensional raster with every "pixel" being a block of 1 meter in length, width and height.

The only command, which my renderer accepts is this:

world.setBlock(x, y, z, type),

where

x,y,z are coordinates and
type is the material of the block (e. g. cobblestone, grass or bricks).

This command puts one block at exactly those coordinates. The space around it remains free. When I want to draw lines, I need to do several calls like that. When I want to fill an area (e. g. a flat roof), I have to do those calls for all of the internal points of the polygon. My renderer doesn't have a "fill polygon" feature, I have to implement it myself.


texts: Getting elevation profile graphs along parallel interpolate lines using ArcGIS 3D Analyst?

texts:
I'm new to ArcGIS. I have used "Interpolate Line" tool in the 3D Analyst Toolbox to get an elevation profile graph.

How can I create several parallel interpolate lines, maybe with specific distance interval, and get their elevation profile graphs?
How can I get the elevation profile graphs from different DEM at exactly the same location?

There's no way that I can change the layer and redraw the same line manually.


texts: How to show the labels of aeroways like Openstreetmap in Mapbox Studio

texts:
I'm using Mapbox Studio to create a Mapbox style.

However, I don't know how to show the labels of aeroways like Openstreetmap.

You can see the labels are showing on the right side.


texts: Filter of pythonaddins.OpenDialog window

texts:
I need to make a pythonaddins.OpenDialog to select a "Folder location"
What is the filternto written in this statement

Folder_Location = pythonaddins.OpenDialog("Please select a folder to store the project and files in", False, "","Select a Folder", ????????,???????)


texts: How do I get the x- and y-coordinate values of points in UTM using QGIS?

texts:
To get to know the size of the study area and how far objects are separated, I need to know x- and y-coordinate values.

Although the $x and $y of Field Calculator returned latitudes and longitudes, x- and y-coordinate values were not obtained when it was in Universal Transverse Mercator (UTM). I think it should be possible to get x- and y-coordinate values since more complicated calculations, i.e. latitudes and longitudes, are possible.


texts: Projecting GCP points into translated file

texts: I have file A.png and file B.png. I translate file A.png so i could translate it into B.png (creating multilayer file). Is there any easy way to get new GCP points coordinates so i could combine those 2 images correctly? I've done so far using only GDAL binaries. Can I somehow use only binaries to do this? If not, is there any easy way do do it?


texts: How to create NSIS installer QGIS with Microsoft Visual C++?

texts:
I want to know, there is a way to generate QGIS NSIS installer with a Release build generated with Microsoft Visual C++?
It is possible to generate the installer in Windows Operating system?


texts: Opening Map Tiles Stored On Google Drive With Google Earth Network Link?

texts:
I'd like to store a superoverlay or tiled map on Google Drive and then open it in Google Earth via a Network Link.

I've been able to store single KML files on Google Drive and then generate a direct download link via a converter and use that successfully in Google Earth as a network link. But I've been unsuccessful in doing this with map tiles/ superoverlays. The network link in Google Maps, which is to the KML, has a green light, but the children files have red lights.

Your thoughts?


texts:
My question is related to this:

https://stackoverflow.com/questions/7927863/location-data-format-for-adehabitat-package

I have a relational database with 'world tiles'. Each tile has a side length of 0.25 degrees (plus top left and lower right lat/long coordinates). This means I have 1036800 tiles (720 for latitudes and 1440 for longitudes). I would like to understand distortion a bit better (sorry about this stupid question). How do you compute the size of each tile? I understand that the distortion is largest at the poles and that the tiles are not necessarily square. Any links/feedback would be very much appreciated. Many thanks!


texts:
I have a simple shapefile (point) that has only ID in the attribute table. I can transform shp to csv but this only gets the ID into the CSV file.

What transformer/settings do I need to read the x,y coordinates of each point in the shapefile and pass them to the CSV writer?


texts:
Have searched the Interwebs for a few hours now but haven't found an answer to a problem which I thought would have a standard solution.

I'm using QGIS to do some analysis on the reach of bus stops in a small town in Sweden with a set buffer of 400m, however I want the results to not have any overlay.
This is the original map

When setting a buffer of 400 meters I get a lot of overlay between stops, however I would like a result similar to this:

Making this "correct" layer took a lot of manual effort, and something I'd rather not do for the remaining 30 bus lines =)

I've tried using Voroni Polygons (with dummy bus stops to expand the reach of the polygons) and a buffer layer, but I see no way to intersect these based on the ID, is there maybe any way to do this?

Or even better, are there any off-the-shelf solutions to creating buffers without overlay which I could use?


texts: How can I translate my Google location history (latitude/longitude) with timestamps to a map that will show me where each pair of coordinates is located?


texts:
I am looking for data showing RARE II Roadless Areas on USFS Lands.

RARE II was an early inventory of Roadless areas that was completed in the late 70s.  This dataset is different than areas that are currently recognized as


texts: 3D features (Polyline, Polygon) intersection library

texts:
For my research I am looking for a (python) library to get the intersection of some 3D features. The features can be Polygons or Polylines. Mostly I need to intersect a plane with my features and get the result.
After some search I came to CGAL, but CGAL build is not successful because of some missing stuff (Eigen) and since I am not very linux, cmake and make literate, I am searching for something else.
(PostGIS  is also OK, but the ST_Intersection seems to work only in 2D).


texts:
If you publish a map with multiple layers using Carto Builder, is there a way for viewers of a map to toggle on/off a specific layer?

For example, imagine a choropleth map with a map marker overlay. While viewing both layers in tandem may be interesting. A user may also want to view just the underlying choropleth map without the overlay.


texts: How to fix Qtcore4.dll can't find procedure entry point?

texts:
Possible Duplicate:
  How to fix QGIS error “Entry Point could not be located”?

I am working with QGIS, python2.5 and pyqt4 for py2.5.
My system variables are as follows:

path:=C:\Program Files\Quantum GIS Wroclaw\apps\qgis\bin;C:\Python25;C:\Program Files\Quantum GIS Wroclaw\apps\qgis;

pythonpath:=C:\Program Files\Quantum GIS Wroclaw\apps\qgis\python;

But when i try to run a .py file through the command shell it gives me this error:

python.exe-entry point not found.
The procedure entry point: ??0QDataStream@@QAE@PAVQByteArray@@H@Z could not be located in the dynamic link library QtCore4.dll

Qtcore4.dll is present in Qgis folder, why does the error occur? and it is not present in system 32 so what do i do then? it was present then i deleted it it still says that.


texts: Which geocoding providers allow storing of geocoded data?

texts:
I'd like to provide a functionality where the user can search for an  address/place on a Leaflet map, and is able to add the resulting marker to the map (which later can be saved on server and downloaded in KML/GPX/GeoJSON format).

The underlying basemap can be anything, it's not guaranteed to be MapBox or Google etc. map.

Now as far as I understand, most geocoding providers prohibit this usage. Google is famous for not allowing it, but even MapBox only allows this on it's enterprise plan, which costs $$$.

Are there any providers which allow you to store the given result from a geocoding request? I believe for places present in OSM, Nominatim and Geonames are a good solutions, but you cannot search for addresses in them.

I'd be especially interested in HERE maps, as they have the highest quality geocoding engine out there (based on this table), but their Terms and Conditions is really not clear about it. It says that layering of their data is not allowed, but it is not clear whether this usage counts as layering.

I'm not looking for a free solutions, but something what is acceptable for a starting project, like Bing, HERE, etc. providing a free usage tier.


texts: Creating Digital Terrain Model from Digital Elevation Model using ArcGIS for Desktop?

texts: I have a DEM and I want to create a DTM from that DEM. Is it possible to do this in ArcGIS?


texts:
I have a value that is passed to a query, and a button that will execute the search function.  The click event doesn't seem to work.  The query works, but it's supposed to run once the button is clicked, and not on map load.

This question Set definition expression from click event  is a bit similiar to my problem, except my query displays.

I was trying to have the user inputted value pass to the query, and run the query when the button is clicked.  I don't want to select features, I only want to show a subset of the layer.

Here is my working code:

......


texts: Unauthorized SQL Query Execution warning- ArcGIS map services

texts:
A security scan of an application I maintain that uses several ArcGIS 10.7.1 map services (hosted with stand alone ArcGIS Server, not Enterprise, on AWS server, with IIS as web server) produced this warning for all of the map services:

Oracle Application Server PL/SQL Unauthorized SQL Query Execution
It is possible to view, modify or delete database entries and tables.

I should note that we do not use Oracle databases. I assume this error is related to ArcGIS Server configuration (Apache Tomcat ?) All the map services are read-only and supported by file gdbs.
The report suggests this to remediate:
Block unauthenticated PUBLIC access to PL/SQL procedures and applications by adding the following rule to the file "$ORACLE_HOME$\Apache\modplsql\cfg\wdbsvr.app":
exclusion_list= account*, sys., dbms_, owa*
Has anyone encountered this problem? I  have never modified any internal configuration to ArcGIS Server.


texts: Filling Colors for Some Countries in Carto

texts:
I imported the "world borders" from the data set library in Carto. I want to fill/highlight some countries (say United States, Africa and Europe) with the same color but not all countries. How can I do this? Is there a way to select and filter the countries?

It is not necessary to do in Carto. If anyone has other tools to do it, I would like to know.


texts:
Is there a simple way to combine KML files into one KML file with Python?

Nothing complex, I just want to append them to each other.

Why I want it to be a script:
Some people in our company create KML files of work site boundaries to display in Google Earth. The idea is to make these all accessible to everyone - my thought is that when a staff member creates a new KML, they can add it to a specific folder on the server. The script will run nightly to combine them, and all staff will have a network link in Google Earth to the folder - thus, they daily get new information in Google Earth.

This was a staff request, so I'm looking into it.

I installed the pyKML library with pip (latest version), but when I try to import it (in Atom with the Script package), I get the following error: ModuleNotFoundError: No module named 'pyKML'

Path appears to be correct.

Can anyone provide assistance with pyKML or an alternate solution?


texts:
I've installed an OpenGTS server and successfully linked it with the GPS2OpenGTS_Pro application for android.

The buses I'm tracking offer 3G with a wireless AP, however whenever the bus stops, the AP is shut down to save battery. The android device has trouble reconnecting all the time, so I've opted for a Raspberry Pi that can connect to the AP through ethernet.

I've got gpsd to read the NMEA output from a BU-353, but I'm stuck trying to figure out how to stream the GPRMC information to my OpenGTS server. I haven't been able to find an application that uploads GPRMC data, and couldn't find a tutorial for how to write a script to do the same (possible using "cat /dev/ttyUSB0" and filtering for just GRRMC lines?). Is there an easy way to do this in python or PHP?

Example HTTP POST:

  http://
  192.3.203.99:8080/gprmc/Data?acct=ACCOUNTNAME&dev=VEHICLENAME&gprmc=$GPRMC,020114,A,3128.7540,N,14257.6714,W,000.0,000.0,131113,,*7


texts: Can this compilation and display of species frequency data be improved? QGIS 1.8.0

texts:
As a newbie on QGIS, I would like to find out if the way I put a species frequency map together could be improved upon.

I had to map the presence or absence of 7 plant species surveyed in 30 2m x 2m quadrats. This data was saved in separate csv files for each species. I then table joined a single csv file to a single vector layer consisting of a 2m grid polygon. The unique ID for the table join involved zooming in on the 2m grid to identify where the survey quadrat had been placed and noting down the polygon identifier number using the identify feature tool. This number was then input as an extra field in the plant species spreadsheet before saving as a csv file.

The final product consists of an background Ordnance Survey 50k Raster map, 7 x 2m grid polygon vector layers each joined to one of 7 x csv files.

Although this has achieved the desired product, I am convinced there must be a slicker way of reaching the same end result, particularly in reducing the number of polygon layers used. Can anyone make any suggestions? Thank you


texts: Iterate incidents for closest facility analysis (Network Analyst) using ModelBuilder?

texts: I have created a model for closest facility analysis. Now I would like to add an iterator to it. Specifically, I need to iterate through 448 "Incidents" (the closest facility model is structured so that it solves from incidents to facilities). Each of the incidents is stored as a shapefile in a folder (one feature per shapefile). So, I would use the "Files" iterator, presumably. But not sure how. Below is a screenshot of the model as currently constructed.


texts:
I'm working with QGIS and I want to reclassify a raster in different ranges.

Using the following statement does not work even recognized as valid

(2950 <= "dem_tallat @ 1"> = 3050) = 3000

Where I can find a document that shows the syntax to use for Raster Calculator?


texts: How do I combine missing data in one table with another?

texts:
I have two tables. I need to join data from table A to table B. Both tables have the same field names and same set of data, except Table A has quite a few rows of data that table B does not have.

The common field between the two is called APN. How can I combine the missing data from table A to table B? I do not want to move over data already inside table B.

ArcGIS software*


texts: Is that possible to use map object ( PHP_MapScript ) instead of map file when using the MapServer?

texts:
In my previous code, I am using OpenLayers WMS to load a raster data through MapServer. As you can see, the mapserv.exe calls a local map file to load the raster data. I just heard PHP_MapScript, since it can convert a map file into a map object.

Is that possible to use this converted map object directly in OpenLayers WMS? Since the object can be easily modified such as "styling".

localWMS = new OpenLayers.Layer.WMS( "Local WMS layers",
          "http://127.0.0.1/cgi-bin/mapserv.exe?MAP=C:/Users/mapserver.map",
          {layers: 'polygon',transparent: true}, {isBaseLayer: false} );
map.addLayer(localWMS);

Someone who is familiar with PHP_MapScript, I am willing to hear from you.


texts:
I need to catch baselayer change events in OpenLayers 3.
In OL 2 I used "changebaselayer" event for map object.

What is a solution in OL3 ?


texts: Using GDAL in python to read elevation at a point from dted files

texts:
How do I get the elevation at a certain point using gdal to read in a dted file? I don't see much useful information out on the web.

I'm sorry I should have been more specific. How would I accomplish this in python?


texts: Adding a graphic line element to open mxd using python comtypes to access ArcObjects

texts:
I attended a pretty awesome demonstration from Jeffrey Barrette at the Esri UC a few weeks ago on creating dynamic tables in ArcMap.  This functionality will greatly benefit the aviation department at my firm.  I have got the scripting all taken care of for creating the dynamic tables in the Layout (using graphic and text elements).  However, this works by cloning a horizontal line, vertical line, and two text elements.

I have already got a working function to add a text element to a map layout via comtypes and ArcObjects, but I cannot seem to get the graphic line part correct. I am using some of the basic functions from a modified version of Snippets (from https://bitbucket.org/maphew/canvec/src/eaf2678de06f/Canvec/Scripts/parco.py).  Here is my function for adding a horizontal line:

def add_horizontal_line():
    GetDesktopModules()
    import comtypes.gen.esriFramework as esriFramework
    import comtypes.gen.esriArcMapUI as esriArcMapUI
    import comtypes.gen.esriSystem as esriSystem
    import comtypes.gen.esriGeometry as esriGeometry
    import comtypes.gen.esriCarto as esriCarto
    import comtypes.gen.esriDisplay as esriDisplay
    import comtypes.gen.stdole as stdole

    # set mxd
    pApp = GetApp()  #current open mxd
    pDoc = pApp.Document
    pMxDoc = CType(pDoc, esriArcMapUI.IMxDocument)
    pMap = pMxDoc.FocusMap
    pAV = CType(pMap, esriCarto.IActiveView)
    pSD = pAV.ScreenDisplay

    # set coords for elment
    pFact = CType(pApp, esriFramework.IObjectFactory)
    pEnv = pAV.Extent
    x = (pEnv.XMin + pEnv.XMax) / 2
    y = (pEnv.YMin + pEnv.YMax) / 2

    # from point
    pUnk_pt1 = pFact.Create(CLSID(esriGeometry.Point))
    pPt = CType(pUnk_pt1, esriGeometry.IPoint)
    pPt.PutCoords(x, y)

    # to point
    pUnk_pt2 = pFact.Create(CLSID(esriGeometry.Point))
    pPt2 = CType(pUnk_pt2, esriGeometry.IPoint)
    pPt2.PutCoords(x + 500, y) # next point 500 units east to form straight line

    # line (from point - to point)
    pUnk_line = pFact.Create(CLSID(esriGeometry.Line))
    pLg = CType(pUnk_line, esriGeometry.ILine2)
    pLg.PutCoords(pPt, pPt2)

    # preset color according to RGB values
    pUnk_color = pFact.Create(CLSID(esriDisplay.RgbColor))
    pColor = CType(pUnk_color, esriDisplay.IRgbColor)
    pColor.Red, pColor.Green, pColor.Blue = (0,0,0) #black line

    # set line properties
    pUnk_line = pFact.Create(CLSID(esriDisplay.SimpleLineSymbol))
    pLineSymbol = CType(pUnk_line, esriDisplay.ISimpleLineSymbol)
    pLineSymbol.Color = pColor

    # create the actual element
    pUnk_elm = pFact.Create(CLSID(esriCarto.LineElement))
    pLineElement = CType(pUnk_elm, esriCarto.ILineElement)
    pLineElement.Symbol = pLineSymbol
    pElement = CType(pLineElement, esriCarto.IElement)
    pElement.Name = 'horizLine'
    pElement.Geometry = pLg

    # add to map
    pGC = CType(pMap, esriCarto.IGraphicsContainer)
    pGC.AddElement(pElement, 0)
    pGCSel = CType(pMap, esriCarto.IGraphicsContainerSelect)
    pGCSel.SelectElement(pElement)
    iOpt = esriCarto.esriViewGraphics + \
    esriCarto.esriViewGraphicSelection
    pAV.PartialRefresh(iOpt, None, None)
    return

I am getting an Error (parameter is incorrect) in this line:

pElement.Geometry = pLg

However, from the API reference it says the ILine2 interface PutCoords method can be used to set coordinates of the line using a FromPoint and ToPoint (IPoint objects).  My understanding of ArcOjbects is not very strong, but it seemed from help docs that this would work.  I even tried to explicitly set the pLg.FromPoint and ToPoint passing in the appropriate coordinates and that did not work either.  Does anyone know what I'm doing wrong here?


texts: Calculating average raster grid values perpendicular to a barrier and transferring values to adjacent cells?

texts:
I was hoping I could get a little help for 2 questions dealing with raster datasets. I have access to ArcMap 10.2 and ArcView 3 to do this. Here is what I am trying to do:

1) I need to calculate an average value for a subset of cells in a line for every cell along a coastline. More specifically, I have a raster grid that lies over the water along a coastline. I need to get the average value of the three cells that lead offshore perpendicular to the coast, this needs to be done for each cell that touches the coast.

2) I then need to move the values calculated above from over water to adjacent cells over land. The over water layer with the values does not extend over the land where I need to transfer the values so I cannot do a simple overlay.


texts:
I wanted to add a GML layer to my QGIS file, but filtered by geometry type. (the file has multiple geometry types in it, with the same layer name, and same layerID.
I came across this post: iface.addVectorLayer on a gml file causes "Select vector layers to add" dialog (pyqgis).
Here they filter it by "|layername="+name
But how can I filter it by geometry type?


texts:
My ultimate goal is to fetch all points from a 4 million record table that fit within a given rectangle and group the results in clusters. All of the queries I have tried up until now take too long. They take upwards of 15 seconds. We are shooting for a few hundred milliseconds max.

I obtained the fastest result from the following approach...

Created a table of 4 million records. Each record contains a geometry point called "location_point". The geometry points are based on SRID 4326.
Created a second column named "snapped_geometry" that is the location_point snapped to a postGIS grid.

set snapped_geometry = ST_SnapToGrid(location_point, 0.2);
Indexed this snapped_geometry column using GIST

create index on tablename using GIST (snapped_geometry)
Ran this query...

explain (analyze) select  count(snapped_geometry) as count, snapped_geometry from contacts_80 where st_contains(st_MakeEnvelope(-95, 30.5, -80, 45, 4326), snapped_geometry)
group by snapped_geometry
This is the response...

Some things I learned from researching the terms in this explain response...
1. The sort information is related to the "group by" clause.
2. The heap scan is related to the "where" clause.
3. Limited work_mem is not the reason for our query taking so long. Initially, there was not enough work_mem to execute the sort in working memory. As a result, the sort spilled to disk. See here. I increased work_mem with set work_mem = '800MB'. This fixed the issue as confirmed by the line "Sort Method: quicksort" in the Explain response.
4. The bitmap heap scan was not lossy. We were initially concerned that our query was lossy because a row in the Explain response displays a "recheck condition". I later learned that this line is in all explain responses even when the bitmap heap scan does not need to recheck the index conditions (i.e. even when the scan is not lossy). The scan only rechecks the condition when the scan is lossy. The absence of the word "lossy" in line 9 of the explain response indicates the scan was not lossy. See here.

Which leaves me still curious how I can speed up this query.
Am I using ST_SnapToGrid incorrectly?
Is there an error in how I created and used the GIST index?
Is it impossible to speed up this query?

I have also experimented with PostGIS's kmeans, clusterDBScan and clusterWithin with no speed advantages.

Other links I used to learn about heap scans and sorting methods...

bitmap heap scans

Sorting methods

Link mentioning count() cannot be used with an index


texts: Snapping lines to points according to attribute field value with PyQGIS

texts:
I would like to snap line features (.shp) to the corresponding point features based on a common attribute field they have (and not according to any distance). I have around 150 different attribute values, and for each I do have only one corresponding point feature but several lines to snap to it.

I thought of doing the following steps to achieve this goal:
1) Creating a list of values containing uniques values from the line attribute field
2) Iterating over the point and the line layer to select features corresponding to each unique value
3) Snap the selected lines to the selected point

I did the first step but I do not manage to code the following steps in PyQGIS. Does someone have any ideas how it could be done?

Here is my network in the first place :

Here the result after snapping the lines to the plot outlet :


texts: Reassigning Raster values to new values in QGIS?

texts:
For example, if the Raster has integer values ranging from -1 to 3 and the value that I am only interested in is 3, I then want to assign values -1 to 2 to 0 and then have the 3 value set to the value of 1. (Similar idea to IDRISI Raster REASSIGN tool.)

I am ultimately trying to find the frequency of a certain cell value in a time series of Rasters of the same extent. To find the frequency of 'cloud' pixels in a certain area over the course of 15 years.

For example, using Raster Calculator to add up the above output Rasters into a single Raster to find the frequency of that cell value since it would have a value of 1 in each Raster.

Additionally, if there is some kind of tool in QGIS that would do a time series frequency comparison on a per pixel basis like the above described that would be good, but so far I have not found such a tool in QGIS.


texts:
I'm a GIS software developer, I've been working predominantly with MapInfo over the last year or so but have now been tasked with assessing QGIS to see if we can perform the same tasks we've developed in MapInfo with QGIS. Or rather compare the performance of the 2. So I'm in the early stages of getting to know how QGIS works, and since I'm used to workspaces in MapInfo I'm trying to understand more about QGIS project files (.qgs).

I've had a scout around but I haven't been able to find any official .qgs documentation, or anything particularly useful elsewhere...does anyone know if this exists?


texts:
I am only looking for a confirmation of what I believe to be true.  If my explaination or conclusion is not accurate I would appreciate it if the reason can be explained to me.

Background Information:

I am attempting to research a method to calculate the Latitude, Longitude, and Height of a point on the planet earth given Range, Azimuth, and Elevation of a point relative to another point where the Lat, Long, and Altitude is known.  I am going to be using a geophysical library that allows me to provide it a geocentric value and it will return a Latitude and Longitude and the Geodetic Datum.

Is WGS-84 a valid datum for a Geocentric datum?
If I provide it a WGS-84 (X,Y,Z) point should it return a WGS-84 Latitude and Longitude?


texts: Reference and Beacon Stations in Vietnam

texts:
I'm not sure if this question can be asked in this forum, but here it goes.
I am on my way to conduct some remote sensing research in Vietnam, which includes logging of coordinates. This should preferably be done with the most exact measure i can get, and I have three hardware options to do this. I have the possibility of bringing a Trimble R8, a Trimble Pathfinder XRS and an ordinary GPS handheld receiver. The R8 works with GPS satellites and reference stations on ground, and the Pathfinder XRS supposedly does the same but with beacon stations instead of reference stations. The handheld GPS works with satellite only and therefore delivers the least precise measurement.

My questions is whether or not reference stations or beacons are located in close vicinity of northern Vietnam, Hanoi. I've looked upon Trimbles website and found a list, but it seems to me that it is deficient in some way. Have you had any experience working with this kind of equipment in Vietnam?

Andreas


texts: How to properly address autocorrelation for logistic regression of spatial data

texts:
Background for the Model:

https://stats.stackexchange.com/questions/268201/how-to-increase-mcfadden-score-and-choose-proper-predictor-for-binary-logistic-r

Data I have:

Forest loss points
Forest boundary
Randomly selected non loss points

(1 & 3 are response variable)

What I need:
I need to minimize the effect of spatial autocorrelation in the regression. For this purpose I need to take points while there is least autocorrelation among the taken points, I think.For this purpose I need to take grids, I think, represent autocorrelation among grid-inside points and least among the points of different grids.So now I need creating those grids for which I need to know the dimension of the grids.

How can I create those grids.

What I have tried so far:
I generated square grids at different spacings and spatially aggregated points by those grids. Finally I ran incremental spatial autocorrelation among those grids at different increment steps but this does not give me a peak point what I would take as grid dimension. Another problem this process creates it gives me a peak point at ~25000m but if I take this distance as the grid dimension then the grid number will be too small and thereby the sampled points would be also too small for fitting the logistic regression.


texts:
When I want to make png tiles that look like ones in osm tile server what I need to do?

Do I need to import data to postGIS db somehow WITH osm Carto style or I can just import it "normally" w/o style
I need to convert CSS style to xml? (If the answer is yes- what is the best way to do this (I've tried with TileMill import))
Last step is to do the "generate_tiles.py" script

I did mange to render files using OSM Bright, but it is not as good as the original osm look for my purpose.

thanks in advance


texts: Cannot find my schemas in QGIS's "Add PostGIS layers"

texts:
When I connect to my PostGIS database using QGIS's Add PostGIS Layers, I see only

public
topology
pg_catalog
information_schema

My other schemas that I have created are not shown here. Anyone else experiencing the same problem? I'm using QGIS v1.8.0-Lisboa, PostgreSQL-9.1 and PostGIS-2.0.1

How should I go about troubleshooting this?


texts: Calculating travel time from residential points to stores via public transit using ArcGIS for Desktop?

texts:
I don't know how to summarize it in just one title, but here is the problem:

My goal is to calculate the travel time from residential points to stores via public transit.  Basically this period of time consists three parts:

The walking time from a residential point to the nearest bus stop "t1", which can be calculated by the Nearest Facility tool from Network Analyst.  The resulting table is called "Resident_Stop", of which the number of records is m;
The bus-riding time "t2" from this bus stop to the nearest bus stop of the target store via bus route, which is a record from Network Analyst's OD Time Matrix table, namely, "Stop_Stop".  The origins and destinations are all bus stops on one single bus route so that the number of records is s*s;
The walking time from the ending bus stop to the store "t3", which is basically the same with the first step.  The table is "Stop_Store", of which the number of records is n.

So the total travel time "t" equals to (t1+t2+t3) and here is where problem occurs: I haven't figured out a way to attach the three tables reasonably to add the times on using a join.  Currently I have an OD Time Matrix table of the walking time directly from residential points to stores via road network (no transit involved) and its number of records is m*n.  I also separatly joined the table "Resident_Stop" in step 1 and the table "Stop_Store" in step 3 to the m*n-sized table based on residential point IDs and store IDs.  So the remaining thing now is to assign the bus-riding time in table 2 based on table 1's stop IDs and table 3's stop IDs to each of the m*n records.  Then I can get the total time.

Should you get any ideas, I'm willing to discuss with you.  I'm now thinking about using Python to handle this but I'm a beginner on that.  If you have any books or tutorials to recommand please tell me.


texts: Problem of downloading Planet data using API

texts:
I am following the method as shown in the Python notebook here to download Planet satellite data from their API.
But, I am getting an error after I run the following piece of code:
import os
import json
import requests
from requests.auth import HTTPBasicAuth

# API Key stored as an env variable
PLANET_API_KEY = os.getenv('PL_API_KEY')

item_type = "PSScene4Band"

# API request object
search_request = {
  "interval": "day",
  "item_types": [item_type],
  "filter": combined_filter
}

# fire off the POST request
search_result = \
  requests.post(
    'https://api.planet.com/data/v1/quick-search',
    auth=HTTPBasicAuth(PLANET_API_KEY, ''),
    json=search_request)

print(json.dumps(search_result.json(), indent=1))

Error:
{
 "message": "Please enter your API key.",
 "errors": []
}

I checked for the API key by doing echo $PL_API_KEY and it produced the correct API key.
How can I figure this out so that I can download data from Planet API?


texts:
SQL within that zone is a point

SELECT
    "public".zona.zts
FROM
    "public".zona
WHERE
    ST_Within (
        ST_Transform (
            ST_GeomFromText (
                'POINT(-30.033514 -51.226828)',
                4326
            ),
            3857
        ),
        zona.geom
    );

I tested in a single step and in separate steps. But the result is always no record.
This pair of coordinates of the point I'm getting via google maps api in a single page geocode addresses, this is my geocode:

<!DOCTYPE html>

    <meta name="viewport" content="initial-scale=1.0, user-scalable=no">
    <meta charset="utf-8">
    <title>Google Maps JavaScript API v3 Example: Geocoding Simple</title>
    <link href="https://developers.google.com/maps/documentation/javascript/examples/default.css" rel="stylesheet">

      <input id="address" type="textbox" style="width:60%" value="Rua Joao Neves da Fontoura,7-Porto Alegre">
      <input type="button" value="Geocode" onclick="codeAddress()">
      <input type="text" id="lat"/>
      <input type="text" id="lng"/>

She gives me coordinates in decimal degrees and my postgis table is in UTM EPSG: 3857 how can I do this conversion? The way I did in the above sql does not work.
Out of curiosity I have sinned a point of a polygon that forms a zone and put as parâmentro and it worked.


texts: CartoDB polygon data for postal codes in Spain

texts:
I have been testing CartoDB for a little while and I can't get maps by postal codes in Spain. The system says "no polygon data available for your selection".

But then I see some maps for Spain including postal code polygons. For example https://dgonzaco.cartodb.com/tables/sppodemos/public/table

Entering the table I see getJSON data containing the polygon data. Is this getJSON data generated by CartoDB, or is it needed to add it manually?


texts: Quantifying mangrove condition over an 8 year period

texts:
Below is the script and I've also put an image of what it's meant to look like and what my image ends up looking like
// The Landsat 8 top of atmosphere reflectance data
var l8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA");

// ROI

// filter the data collection
var filtered = l8.filterDate('2013-04-01', Date.now())
        .filterBounds(roi);
print(filtered);

// Get a cloud free image
var cloudFreeImage = filtered.sort('CLOUD_COVER').first();

// add the cloud free true colour image to mapping layer
var rgb_vis = {min: 0, max: 0.3, bands:['B4', 'B3', 'B2']};
Map.addLayer(cloudFreeImage, rgb_vis, 'True color');

// add the cloud-free false-colour image to the mapping layer
Map.addLayer(cloudFreeImage, {min: 0, max: 0.3, bands:['B5', 'B4', 'B3']}, 'False Color');

// add the median image to the mapping layer
Map.addLayer(filtered.median(), rgb_vis, 'median RGB');

// compute NDVI and add the NDVI image to the mapping area
var ndvi = cloudFreeImage.normalizedDifference(['B5', 'B4']);
Map.addLayer(ndvi, {min: 0, max: 1, palette:['brown','yellow','green']}, 'NDVI');

// add a new ndvi band a single l8 image
function addNDVIband(anImage) {
  var ndviTemp = anImage.normalizedDifference(['B5', 'B4']);
  return anImage.addBands(ndviTemp.rename('NDVI'));}

// test the function on a single image
var ndviAdded = addNDVIband(cloudFreeImage);
print ('Before adding NDVI', cloudFreeImage);
print ('After adding NDVI', ndviAdded);

// map the function to the image collection
var ndviMapped = filtered.map(addNDVIband);
print ('After mapping NDVI', ndviMapped);

// apply temporal reducer to the NDVI band that has been added to the image collection.
Map.addLayer(ndviMapped.median(), {bands: 'NDVI', min: 0, max: 1, palette:['brown','yellow','green']}, 'median ndvi');

/*// chart the ndvi over time
var ndviChart = ui.Chart.image.series({
    imageCollection:ndviMapped.select("NDVI"),
  region: roi,
  reducer: ee.Reducer.median(),
  scale: 90,
  xProperty:'system:time_start'});
print(ndviChart);
*/

[][1
<img src="https://i.stack.imgur.com/r7Fp9.png" alt="what mine looks like with the code I'm using
" />


texts: Looking for examples on how GIS techniques facilitate utility surveying and management

texts: Anyone get some ideas about the interface between underground utility and GIS? Utility surveying is quite a new topic to me and I want some examples on how GIS techniques facilitate utility surveying and management.


texts:
I am trying to run the following PostGIS query using CartoDB's SQL API:

SELECT
the_geom,
ST_Area(ST_Intersection(the_geom::geography,ST_Buffer(ST_SetSRID(ST_MakePoint(-52.4379, -12.6235),4326)::geography, 25000))) / 1000000 AS polyarea FROM
    dataset1
WHERE
    ST_Intersects(the_geom::geography,ST_Buffer(ST_SetSRID(ST_MakePoint(-52.4379, -12.6235),4326)::geography, 25000))

The same query works with two other datasets (one of which has 735,000 rows), but returns a

  "row too large error"

for this particular dataset (which "only" has 160,400 rows).

Update: Is it possible that the features are too complex for the calculations?


texts: Exporting a QGIS layout with Addition on layers

texts:
I've been playing around with an approach in QGIS that uses Addition between the layers I'm mapping to build a tri-variate composite colour map - an example is here in the QGIS Flickr pool:

http://www.flickr.com/photos/97267218@N06/12626847074/

To achieve an output I thought I had to set the background of the image to black - in order to draw out the colours. However I've been asked (demanded!!!) to set the background to white.

Thinking about it a little, one approach that would seem to work would be to set the background of the Composer to Alpha 0.

This is a screen grab of an map produced in a similar manner. In the composer I set the background with an Alpha of 0 - meaning that it doesn't influence the layers. Great!

However when I export the final image in any format (including Pdf and SVG) I only get a white page!

There has been a similar question asked about this in ArcGIS (Exporting Graphic with transparency in layout to PDF using ArcGIS for Desktop?) - and the suggestion there was to change the graphic to a feature and insert that .... something I don't think is possible with QGIS ...

So - has anyone either any suggestions how this might be overcome? How can I export an image which is using Addition and Alpha without needing a black background?

More fundamentally ... can anyone explain (in relative layman's terms!) why the export isn't working in the first place?


texts: How to do a spatial join in PostGIS that returns on the minimum value of point intersecting polygon?

texts:
I have a polygon table (buffer) and a point table (dem_points), and I want to do a spatial join on them.  I need to join the point with the minimum elevation (dn) to its spatially corresponding polygon in PostGIS.

In the above example I want my output to be the polygon with osm_id of'1234' joined with the point where dn=4200.
I have this query but something is off:

SELECT osm_id, gid, MIN(dn)
FROM (
      SELECT a.id AS osm_id, b.id AS gid, b.dn AS dn
      FROM
         buffer AS a, dem_points AS b
      WHERE ST_Within(buffer.geom, dem_points.geom)
     )
GROUP BY (osm_id, gid);`

With this query I receive this error:

ERROR:  subquery in FROM must have an alias
LINE 2:      FROM (
              ^
HINT:  For example, FROM (SELECT ...) [AS] foo.
ERROR: subquery in FROM must have an alias
SQL state: 42601
Hint: For example, FROM (SELECT ...) [AS] foo.
Character: 40

How do I need to structure my query to spatially join the point with the minimum elevation to the polygon it falls within?


texts:
I am using Leaflet 1.3.1, for show a map in a blog made with Laravel and Vue.js

The screen it is showing in grey, however, the markers and controls are showed.

When I check Inspect Elements in my browser, and check the layer , I can see that the .png of the map are loader, if I push it, the Tile or PNG is showing correctly and I can see this part of the map.

This is my code:

layout/app.blade.php

  ..../
<link rel="stylesheet" href="https://unpkg.com/leaflet@1.3.1/dist/leaflet.css" />

/.....

Home.blade.php

  extends('layouts.app')
  @section('content')

@foreach($articles as $article)

@endforeach
@endsection

I could check in this forum and checked answer like this:

Leaflet output is grey

However, if you check my code, you can see that I tried to resolve it, but it not working yet.


texts: Retrying failed TileWMS tile load in OpenLayers

texts:
I have a map in OpenLayers with TileWMS source that's not very reliable. At almost every full page load there are a few tiles that fail to load with status 503 - Service Unavailable. When I try to load failed tile a bit later 'manually', it loads without a problem. That lead me to the idea of reloading failed tile if first try is not successful.
I couldn't find any OpenLayers 'official' way of doing that, so I tried to achieve that using tileLoadFunction, where I initiate reload of the tile with 2 seconds delay and before reload check if tile was already loaded successfully. Reload is triggered by first assigning dummy image source for tile and then the original one. Successful loading of the tile is checked by hooking on tileloadend event and setting custom flag evt.tile._loaded = true.
Here is the relevant code:
proj4.defs('EPSG:3912',
           '+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=500000 +y_0=-5000000 +ellps=bessel +towgs84=682,-203,480,0,0,0,0 +units=m +no_defs');
var proj3912 = ol.proj.get('EPSG:3912');
proj3912.setExtent([374371.84, 30513.32, 624119.18, 195517.48]);

var tileGridResolutions = [1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1, 0.5];

var dofLayer = new ol.layer.Tile({
  source: new ol.source.TileWMS({
    url: 'https://prostor4.gov.si/ows2-m-pub/wms',
    attributions: ['GURS'],
    params: {
      'LAYERS': 'SI.GURS.ZPDZ:DOF050_D48',
      'FORMAT': 'image/png',
      'VERSION': '1.3.0'
    },
    tileGrid: new ol.tilegrid.TileGrid({
      extent: proj3912.getExtent(),
      resolutions: tileGridResolutions,
      tileSize: 256
    }),
    tileLoadFunction: tileLoadFunction,
    projection: proj3912
  }),
  title: 'Ortofoto',
  type: 'base'
});

function tileLoadFunction(imageTile, src) {
  imageTile.getImage().src = src;
  if (typeof imageTile._retry == 'undefined') {
    imageTile._retry = 0;
    imageTile._loaded = false;
  }
  if (imageTile._retry < 3) {
    imageTile._retry++;
    setTimeout(function() {
      if (!imageTile._loaded) {
        imageTile.getImage().src = 'img/empty.png';
        tileLoadFunction(imageTile, src);
      }
    }, 2000);
  }
}

dofLayer.getSource().on('tileloadend', function(evt) {
  evt.tile._loaded = true;
});

Unluckily this method works only in old venerable IE11, but not in Firefox and Chrome. There is no error, it just seems like image is not transferred from image element to canvas. Since I work in ES5 JS, I cannot debug OpenLayers internals to see where the problem is.
Is there any way to make the above method work in Firefox and Chrome or some other approach to try reloading failed TileWMS tile?


texts:
I have two layers, the first active in the atlas and the second in the table of contents. Both are related with a id "Refcat".

I want show only those actives in the atlas with their own symbology.

The first

$id = @atlas_featureid

And the second is coverage.... How to do this?


texts: Scripted generating then loading rasters to PostGIS - raster2pgsql phantom

texts:
(Ubuntu server 18)
It seems that raster2pgsql is included in the postgis package but I can't find it separately to download either source or binary. Supposedly there's an executable somewhere in the postgresql lib folders. After reading some things I have the idea that it's not default in the installation or when compiling and maybe I won't find it if I used the extension script distributed with postgresql. How do I add it?
That said, the intention is to generate raster data and load it using Python scripts, so a gdal driver might be more convenient since we are already using gdal.
The script runs on a separate machine from the pg instance, if that matters.
What's the most straightforward way to pipe large volumes of raster data into postgis on a remote box?


texts:
I'm currently working on displaying twitter geo data (currently just long/lat -> cartesian without any analysis).

This is the world map (source: naturalearthdata): http://dl.dropbox.com/u/14747311/world.png

This is the rendered data: http://dl.dropbox.com/u/14747311/overlay.png. The mapping is done via mercator projection (with R = 1):

x = ((long + 180)/360)*x_dim
y = (y_dim/2) - ((arsinh(tan(lat))*(y_dim/2) (lat converted to radiant of course)

Now, the shape is apparent, but of course, the projections don't match. How can I fix that? (programming experience is plentiful, so it's more the theoretical side I'm concerned with).


texts: Interpolating Maps - Statistical Learning Techniques vs Spatial Statistics Techniques

texts:
I have recently begun experimenting with the interpolation of variales over a whole country.

I have good sized training set where I have several dependant variables (continuous and categorical) associated to many independant variables. The independant variables are available for the whole country.

I'm not at all an expert in spatial statistics so my first approach was to train machine learning models like random forests, svm, gradient boosting, etc Which are behaving very well.

But I have been finding out that there are default methods for interpolating available in gis packages. Kriging. I played a bit with Kriging was not satisfied. Could co-kriging - regression-kriging get to better perfomances than these before mentioned methods (rf, svm)?

Also, I have noticed that there are methods which seem to be tailored for these problems like geographically weighted regression. I am just wondering if my approach is somewhat off if there are methods specifically designed for these tasks.

Thank you.


texts: How to apply schema changes on Esri geodatabase two-way replica set?

texts:
I have an Esri geodatabase two-way replica set. I need to update the schema on both databases in the replica set. Esri describes doing this using ArcGIS schema tools (compare/export/import replica schema - see here).

I have control of both databases in the replica set and I would prefer to make a python script, which updates the schema and run that on both databases in the replica set.

Are there potential problems with applying the schema changes to both databases instead using Esri's schema tools?


texts:
Does anybody know a trick to create an antipodes map with QGIS?
That is a map of the world on one layer, and a map of the world on another layer, showing the antipodes of each point (practically a world map inverted north-south and east-west).
I did not find any shapefile with the inverted world in the net, so I guess I will have to play with the projection, which is a very complicated issue.


texts: Why is my Points Distance result a decimal instead of miles

texts:
I have a bunch of points representing residential addresses of employees and two points representing office buildings.

When I use the points distance function, the distance comes up as a decimal and not as the correct mile distance.

Any ideas why that might be?


texts: How to speed up Postgis raster layer display in QGIS?

texts:
i've been importing some digital orthophoto geotif files into postgis 2.0.1 via raster2pgsql tool (i also made 3 levels of an overview table - pyramid). Then i tried to import the layer into qgis via the postgis load raster plugin. It works, but the speed seems so slow its unusable at the moment.

I found a post somewhere about this being a GDAL issue, but i couldnt find if a solution has been found yet. Did anyone get across this?

Thanks


texts: Is it possible to merge two geoserver layers into one (merging) and export as KML?

texts:
I have used the KML reflector in order to get a geoserver layer into KML.
Now what I want to do is to merge two layers (which might have different geometries) and export them into a KML file.
Is something like this possible with openlayers or geoserver api?


texts: How to dissolve geodatabase map across two subdivision levels in Arcgis 10.2

texts:
The existing polygon database contains attribute table with the following structure, and has three levels; country, state, and city (and many other countries not listed below).  I want to dissolve the map at country level, but for some specific countries. I want further dissolve them at state level.

Thanks.

ID_0    country ID_1    state   ID_2    city    ID_state    ID_city
48      China   01      Anhui   08      Hefei   4801        480108
48      China   01      Anhui   08      Hefei   4801        480108
48      China   01      Anhui   15      Tongling4801        480115
48      China   02      Jiangsu 14      Suzhou  4802        480214
48      China   03      Fujian  15      Fuzhou  4803        480315
.
.


texts:
I have created two labelled vector points on an OpenLayers map. I'm trying to change the label text dynamically. I tried changing the values where the stylemap gets the values from and calling layer.refresh() but it does not do anything. Any ideas?

Here is my code in its entirety. All you have to do is change the OpenLayers JavaScript reference to your local file and it will work in your browser.

    <title>OpenLayers Labeled Features Example</title>
    <link rel="stylesheet" href="../js/OpenLayers/theme/default/style.css" type="text/css" />
    <link rel="stylesheet" href="style.css" type="text/css" />

OpenLayers Refresh Labelled features example

        <input type="button" value="Switch Labels" onclick="Javascript: switchLabels();" />


texts:
I work with 72 MXD's files, and i need to turn off the grid in order to export it without the grid.

I don't want to do it manually, neither using the "Tinytask" software, as mentioned in Changing grid interval using ArcPy?

Is it possible to do this using python stand-alone code?


texts: Comparing two TINs created using ArcGIS for Desktop?

texts:
I have the following two files:

LiDAR las file1
Let's give it a name LiDAR las file2, contains the same file1 data, where number of points randomly removed.

So if file1 contains X number of points, file2 would contain Y number of points, where Y is less than X.

Now, I am using ArcGIS "LAS to Multipoint" to create from both files above their shape counterparts. Next I am creating two empty TIN(triangulated irregular network) files using ArcGIS. Then I edit both TINs, let's say tin1 and tin2, with the appropriate input feature classes, i.e. file1 to tin1 and file2 to tin2.

Is there any tool/package/library or any software that would let me calculate the RMSE(root mean square error) or any difference information that would give me information on the accuracy of the second file, i.e. file2, with respect to file1? I tried to use the TIN difference but it gives system error and cannot obtain the result. I would really appreciate if anyone could show how to obtain accuracy information and interpretation of the results obtained.


texts:
I would like to have a clean viewshed analysis of an openDEM region.
I've got this elevation TIFF:

But after running the Viewshed "Visibility analysis" plugin, I obtain this result:

It's the same with the GRASS viewshed tool.
What produces these NULL pattern values on my result raster? The altitudes are consistant in my openDEM raster.


texts: Driver error : DTM to raster using GDAL_GRID

texts:
I am looking to interpolate a DTM (eventually over 4,000 DTMs) that is in a personal geodatabase format.  I would like to do this in as few step as possible, eventually writing a batch script to do them all.

I was under the impression that I could use gdal_grid to accomplish this but I am getting the following error:

ERROR 1: Unable to initialize ODBC connection to DSN for DRIVER=Microsoft Access
Driver (*.mdb);DBQ=17509200.mdb,
[Microsoft][ODBC Driver Manager] Data source name not found and no default drive
r specified
Unable to open input datasource "17509200.mdb".
Unable to initialize ODBC connection to DSN for DRIVER=Microsoft Access Driver (
*.mdb);DBQ=17509200.mdb,
[Microsoft][ODBC Driver Manager] Data source name not found and no default drive
r specified

after using the following command:
    gdal_grid 17509200.mdb test1.tiff -zfield ELEVATION -ot Float64 -l 17509200

I am using GDAL 1.9.2 on a Windows machine.

My version does support reading of "Pgeo" so I thought GDAL would be able to use the .mdb as an input, and output a raster. But apparently I have to specify a driver (I'm not sure why it can't use PGeo?).  I have seen on other forums Linux users using various MDB driver extensions.  Does anything exist similar to that for Windows? Or can I specify the PGeo driver in GDAL as it is now? How do you specify a driver in GDAL?

Any suggestions?


texts: Clip Dataframe Arcgis for data driven map

texts:
I'm trying to do a map of several provinces bordering each other. The map needs to focus on each province and not the neighbors. therefore when using the data driven map function i would like to clip the dataframes extent to match the province in which is being exported...

obviously we can do them individually however one click functionality would be better.


texts:
I believe I understand the shapely methodology fairly well at this point.  However, I'm struggling with getting my holes to populate to the polygons.  I've tried a number of approaches.

What I'm doing is taking a buffer feature that I then want to simplify the edges.  So, I take the exterior and the individual interior LineStrings and split them by the points that were used to create the buffer.

All of this works fine.  However, I then want to take that exterior and the separate interiors and create a single Polygon feature with the holes.  Here's some of the various methods that I have tried:

polygonExterior = Polygon(self.SplitLinesByPoints(inputPolygon.exterior,...)) # returns a LinearRing
polygonInteriors = []
for i in range(len(inputPolygon.interiors)):
    polygonInteriors.append(self.SplitLinesByPoints(inputPolygon.interiors[i],...))

#when I check polygonInteriors, it has 3 valid LinearRings

#Trial 1
newPolygon = Polygon(polygonExteriors, polygonInteriors)
# returns a polygon with the correct exterior and no holes (should be 3 in my test dataset)

#Trial 2
interiorCoords = []
for interior in polygonInteriors:
    interiorCoords.append(list(interior.coords))
newPolygon = Polygon(polygonExteriors, interiorCoords)
# same result as Trial 1

#Trial 3
newPolygon = Polygon(polygonExterior)
for interior in polygonInteriors:
    newPolygon = newPolygon.difference(Polygon(interior))
#same result

#Trial 4 based on first received answer
newPolygon = Polygon(polygonExterior, [[pt for pt in inner.coords] for inner in polygonInteriors])
#same result

# testing output
print(str(len(newPolygon.interiors)))
# gives 0

I output the list of list of tuples that are passed in trial 4 and it looks like this:

[[(3149326.399938697, 13547034.20647637), (3149326.3938959567, 13547034.177421806),...,(3149326.399938697, 13547034.20647637)], [(3149339.754347628, 13546944.667812834), (3149339.6827200023, 13546944.533919774),...,(3149339.754347628, 13546944.667812834)], [(3149334.7216034834, 13547168.84576934), (3149334.729913973, 13547168.803259838),..., (3149334.7216034834, 13547168.84576934)]]

I also just wanted to confirm that the interiors were where they should be. So I created shapefiles for each of the interior circles as well.  This is what I expect to get:

As you can see, I have the actual output with the three holes as separate shapefiles.

But, what I get is this:

That's just the TestInterior2... file that should have the holes.

I'm just really at a loss as to what to try next.

Anyone have any ideas what I'm doing wrong or how to get the holes into the Polygon?


texts: Flow Accumulation and Flow Path Error from Flow Direction Raster in ArcMap?

texts:
I am trying to generate a flow accumulation and flow length raster from flow direction.

I have successfully generated a flow direction raster from a DEM that was clipped by the 'Clip (Data management tool)'. This was done using the hydrology toolset in spatial analyst tools, in ArcMap 10.2. See image.

When I attempt to generate the flow accumulation or flow length rasters in ArcMap 10.2 I am receiving the following error message:

  'Output accumulation raster' <empty>.

Has anyone had a similar problem and rectified it? Could it be something to do with the way the data has been clipped? I only want to generate flow accumulation and flow lengths for the clipped area.

This problem was solved by following the advice of FelixIP 53 in comments.


texts: I have encountered a GeoTIFF raster files with EPSG:4326 (WGS84), i.e. no projection defined! Is this normal? I would expect raster to be in Projected Coordinate System. Raster is a grid of pixels so I expected this would be defined in the plane.


texts: How to correct a wrong projection of shapefile?

texts:
I have a shapefile with UTM zone 16. This projection is wrong. It must be UTM zone 32. How can I correct the projection?

When I convert (f-tool) to UTM zone 32, all the vector are rotated. The property of the layer shows funny numbers that is not correct:

xMin,yMin -4417494.80;6867783.98 : xMax,yMax -4416554.49;6869279.82

Correct property should be:

xMin,yMin 3410596.14;5317812.48 : xMax,yMax 3411876.35;5318694.50


texts: QGIS: building pyramids faster in properties compared to gdaladdo

texts:
Does anybody know why building pyramids in QGIS using the "build pyramids" button from the properties menu (tab: pyramids) is a lot faster than command line gdaladdo or Raster->Misc->Build Overviews?

To rule out some questions that may arise:

I didn't build less pyramids with the properties dialogue
I didn't choose different methods


texts: How to change default image in tree node with GeoExt

texts:
I'm using GeoExt to create a tree. In each node I want the corresponding image at left instead of the default image   like this site. How can I do it?

Edit

I found this example but I don't know how can I add childrens

Edit 2

This is part of my code:

var bpn = new OpenLayers.Layer.WMS("bpn",
    url,
    {
        LAYERS: [ 'Año 2004', 'Año 2005', 'Año 2006', 'Año 2007', 'Año 2009'],
        format: "image/png",
        transparent: "true",
        projection: 'EPSG:4326'
     },
     {
         buffer: 0,
         displayOutsideMaxExtent: true,
         isBaseLayer: false,
         displayInLayerSwitcher: false,
         yx: {'EPSG:4326' : true}
      });

var treeConfig = [
    {
        nodeType: "gx_layer",
        layer: bpn,
        isLeaf: false,

        loader: {
            param: "LAYERS",

        },
        expanded: false
    }];

tree = new Ext.tree.TreePanel({
    border: true,
    region: "west",
    title: "Entre Ríos",
    width: 250,
    split: true,
    collapsible: true,
    collapseMode: "mini",
    autoScroll: true,

    loader: new Ext.tree.TreeLoader({
        applyLoader: true,
        uiProviders: {
            "layernodeui": LayerNodeUI
        }
    }),
    root: {
        nodeType: "async",
        children: treeConfig
    },
    rootVisible: false,
    lines: false,
});

Each layer (node)'Año 2004', 'Año 2005', 'Año 2006', 'Año 2007', 'Año 2009' has a different image. What should I change in treeConfig?


texts: How to convert float raster to vector with python GDAL

texts:
I have float raster and now I want to convert it to vector. How is it possible with the Python GDAL library?

I have tried with gdal_polygonize.py of GDAL utilities on the command line and it worked excellently. But this utility is based on GDALPolygonize() of C++ library and I want to C++ method GDALFPolygonize() to be used instead, which manipulates float raster data as far as I know.


texts: Updating Pyproj syntax results in other values


texts: SELECT * will not work with table containing geodata

texts:
The following will not run... I get ORA-00932 error every time..
$sql = 'SELECT * FROM app.table_1';

however, this will work (specifying a column name, i.e. objectid):
$sql = 'SELECT objectid FROM app.table_1';

Seeing as I have a table with ~100 columns I do not want to manually hardcode out a list of the columns for this query.. anybody experience this?
legacy question: https://stackoverflow.com/questions/63562690/cant-select-ora-00932-inconsistent-datatypes-expected-char-got-adt


texts: ArcGIS loads splash screen then quits (possible Python issue)

texts:
I have a student license of ArcGIS Desktop 10.6.1, working on Windows 10.

When I try to open ArcMap, the splash screen loads, but then the program never actually opens. It doesn't generate any crash report, it just closes (have checked many places for error reports).

Other ArcGIS programs work (i.e. ArcCatalog), but close instantly if I open the Python interpreter.

If I install Arc without Python, the program works and the Python terminal will open, but it doesn't actually work. I.e. even the most basic print statement does nothing. Otherwise, I get the splash screen and 'crash'.

I suspect that it has something to do with Arc trying to use the 64-bt python version instead of the 32-bt version, but the instructions I found were unclear and I'm not sure how to fix this.

Things I have tried in Arc:

Repairing Arc; modifying by removing and then adding Python again workaround 2, 3 3
Full uninstall/reinstall of Arc
Modifying registry (python) workaround 1
Modifying registry by deleting software>ESRI 4
Updated graphics card driver (have Intel UHD Graphics 620)
Deleted anything related under AppData\Roaming
Modified/deleted normal.mxt 4
Verified contents of .pth (didn't need any changes) 5
Checked environment variables 5

Other things I tried on my machine:

Verified all computer updates are installed and ran a tune-up
Updated graphics card driver (have Intel UHD Graphics 620) 6
Checked for any driver conflicts, found none 7

My User PATH:

My system PATH:

I also added: C:\Python27\ArcGIS10.6\python.exe

I've tried contacting ESRI support but because it is a student version, any request has to go through my prof (who is being unresponsive). I need both ArcMap and Python/ArcPy to work, and I can't remove Python 3.7 from my machine because I need it.

Update: if I rename my Python37 folder Arc starts working - it opens and the python window works. So clearly there's something directing ArcGIS to 3.7 instead of 2.7. It's not great to have to rename my Python37 directory every time I want to use Arc though - is there another way??


texts:
I'm interested in creating a component / plugin for ArcMap that works something like ArcGIS online. Ideally a small web app would boot up within a window in ArcMap which would be able to perform functions such as adding a layer and other standard tasks.  In the first instance it would be a search tool, but may perform other activities in the future.

Does anyone know if this is possible within ArcObjects?  Is there an API, perhaps with a set of functions that can be called in JavaScript that can do this?

Thanks

Will


texts: QGIS: Symbolize by date field based on current date

texts:
In QGIS 3.2.3: I have a layer/Shapefile with a date field (type DATE) that comes from a joined table. I would like to symbolize the points which are older than 10 days (based on the current date) in another color than points which are younger than 10 days (based on the current date). It would also be possible not to display the younger points.

Is there a way to use the current date to make such visualisations? Maybe with a rule based symbology: I thought in using an expression with the now() function, which returns the current date and time and then remove 10 days from it, but couldn't get it going. Am I missing something?


texts:
As you can see on the screen shot below, I'm trying to create a 3 km buffer around the roads.
After clicking OK, the buffer is successfully created and added to the table of content but nothing shows or changes on the actual map.

Can anyone help me with it?


texts:
Is there any free solution for geocoding UK with Java?

I would like to avoid APIs, like for example Google's Geocode service as I have over 1 million addresses to geocode.


texts: Using Proj.4/GDAL to convert lat,lon to i,j

texts:
I'm working with a data set in GRIB2 format and need to determine the grid cell index (i, j) corresponding to a given lat,lon.
I have the projection definition in the form of projparams as follows:

{'a': 6371229,
 'b': 6371229,
 'lat_0': 25.0,
 'lat_1': 25.0,
 'lat_2': 25.0,
 'lon_0': 265.0,
 'proj': 'lcc'}

which gives me the following srs:+units=m +a=6371229 +b=6371229 +lon_0=265.0 +proj=lcc +lat_2=25.0 +lat_1=25.0 +lat_0=25.0

And I know dx, dy, number of grid points in each direction etc.

So I was thinking of using Proj.4 (because it has JavaScript bindings) (and/or GDAL) to convert it to i,j. I think I have read somewhere that one could do it this way by defining some kind of equidistant metric map projection that would use i, j as x, y coordinates and do a simple coordinate transformation, but I'm struggling with how to actually define such a projection in proj.4 terms.


texts: Publishing a three band ArcGIS ImageService and encoding tiles with LERC

texts:
I'm trying to publish an ImageService for a three-band mosaic dataset and serve tiles using the LERC format. My workflow is as follows:

Create the mosaic dataset: right click on a file geodatabase --> New Mosaic Dataset --> Pixel Properties --> set Number of Bands to 3
Publish as image service: once created, right click on the mosaic dataset --> Share as Image Service --> Caching --> Select "Using tiles from a cache"
Specify tile format: After specifying that this dataset will be cached, the Advanced Settings menu appears below Cached. Within this menu, LERC is not available from the Tile Format dropdown

LERC is available if I follow the above, but with a 1 band mosaic dataset.

Is this option be available in 10.4? Or will it be in 10.5? I'm able to read three band images with LERC encoding using the Export Image endpoint on an ImageService, but need to be able to consume these as tiles as well.


texts: How to fix "Error processing specified tolerance!" in Polygon Centroids tool?

texts:
I am using QGIS 1.8.0 to find the centroids of 92 features (i.e. smaller polygons) in a polygon layer.

However, while using (Vector -> Geometry Tools -> Polygon centroids), I encountered the following error:

"Error processing specified tolerance!
Please choose larger tolerance..."

Can someone please explain to me what went wrong and how I can fix it?

Thanks!


texts:
Is there any real advantage in using L.GeoJson instead of Polyline when drawing paths in a map powered by Leaflet?

I know that GeoJson can attach additional options. It is also possible with Polyline but one has to extend it or have a wrapper object.

Our project needs dynamic addition and removal of paths. Can GeoJson still handle it?


texts: > - operator and "Shape_Area" field in a geodatabase

texts:
this might be a pure python problem or a problem with arcpy. please help me out:

    min_area = float(GKPROZ_functions.min_of_field(out_name,"Shape_Area"))
    print area_max
    if min_area > area_max:
        print "success!"
    else:
        print "failed."

    if min_area > 100:
        print "success!"
    else:
        print "failed."

    print str(min_area)

output:

100
failed.
success!
100.018289077

min_of_field returns the minimum value of the field "Shape_Area" of the whole feature class.
it is (apparently) 100.018289077. however if i check that with area_max (supposedly 100) it says it's smaller (failed!). if i check it with 100, it says it's bigger(success!).
what mistake am i making here? thanks in advance.


texts:
Trying to render a circle and buffer from the WFS link generated by me, with GeoJSON format using Leaflet in an app and QGIS on desktop to render it, the image below was the result.

How do I fix this?


texts: Accessing values in GIS based on condition using NumPy?

texts:
I wrote the following code for GIS analysis, in which I am trying to find the sum of all pavement area values that meet the condition where Application == "Micro Surface" & CouncilDistrict == 1.

I have more than 85,000 records in the attribute table with any of six different values for the field named Application and values ranging from 0 to 10 for CouncilDistrict. I want to find the sum of values for the Pvment_Area field that meet the condition. I keep getting the following error:

  ValueError: The truth value of an array with more than one element is
  ambiguous. Use a.any() or a.all()

How do I extract values using NumPy?

import arcpy
import numpy
input = "C:\Temp\Street1.mdb\Scenario_78M"
arr = arcpy.da.FeatureClassToNumPyArray(input, ('Pvment_Area', 'Application', 'CouncilDistrict'))
x=(arr[arr['Application'] == "Micro Surface" AND ['CouncilDistrict'] == 1]]["Pvment_Area"].sum())
print'Pvment_Area is',x


texts: Handling PostgreSQL IntegrityError from PyQGIS

texts:
I am developing a QGIS plugin, in which I want the user to be able to add/delete features from some layers.
Here is the code I'm using to delete a feature with id = trayId

try:
   self.traysLayer.startEditing()
   self.traysLayer.deleteFeatures([trayId])
   self.traysLayer.commitChanges()
except Exception:
   self.iface.messageBar().pushMessage("Error", "The selected tray can't be deleted because it is in use.")

The problem is sometimes the feature to be deleted is referenced by another table using foreign key, so I'm getting the following PostGIS error:

ERROR: update or delete on table "tray" violates foreign key constraint

I have tried adding try/except to handle all exceptions, but it doesn't succeed to handle that error, I am still getting the same error.

How can I handle such error?


texts:
How to find the Query Builder in the freshly released QGIS version 3.0?

In previous versions (2.x), it could be found as mentioned in Displaying only selected features on map in QGIS?


texts: I have more than 50 raster files. I have used to mosaic to new raster tool on Arcgis and I wonder if there's anyway I can select all rasters instead of having to click one by one to add them up.


texts: OpenLayers call ArcGIS WFS: no feature returned

texts:
I spent a long time looking for an example using OpenLayers to call ArcGIS WFS, however, almost none is found... currently the code I wrote is based on the example of GeoServer, however, no feature is returned...

new OpenLayers.Layer.Vector("ArcGIS WFS layer", {
    strategies: [new OpenLayers.Strategy.BBOX()],
    protocol: new OpenLayers.Protocol.WFS({
        version: "1.1.0",
        srsName: "urn:ogc:def:crs:EPSG:6.9:26912",
            url: "http://XXXX/Geodatabase_SDE/geodataserver/WFSServer?request=getFeature",
        featureType: "TEST_SAMPLE",
        featureNS: "http://XXXX/Geodatabase_SDE/geodataserver/WFSServer",
        geometryName: "Shape"
    })
});

I am wondering if it's doable in OpenLayers?

UPDATE
Request XML:

<wfs:GetFeature service="WFS" version="1.1.0" xmlns:wfs="http://www.opengis.net/wfs" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
  <wfs:Query srsName="EPSG:26912" typeName="feature:TEST_SAMPLE" xmlns:feature="http://XXXX/Geodatabase_SDE/geodataserver/WFSServer">
    <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
      <ogc:BBOX>
        <ogc:PropertyName>
          Shape
        </ogc:PropertyName>
        <gml:Envelope srsName="EPSG:26912" xmlns:gml="http://www.opengis.net/gml">
          <gml:lowerCorner>
            447517.58546862 4364229.0103902
          </gml:lowerCorner>
          <gml:upperCorner>
            455008.41453138 4384331.9496098
          </gml:upperCorner>
        </gml:Envelope>
      </ogc:BBOX>
    </ogc:Filter>
  </wfs:Query>
</wfs:GetFeature>

Response XML:

<wfs:FeatureCollection xmlns:Geodatabase_SDE="http://XXXXX/Geodatabase_SDE/GeoDataServer/WFSServer" xmlns:gml="http://www.opengis.net/gml" xmlns:wfs="http://www.opengis.net/wfs" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://XXXXX/Geodatabase_SDE/GeoDataServer/WFSServer http://XXXXX/Geodatabase_SDE/geodataserver/WFSServer?request=DescribeFeatureType%26version=1.1.0%26typename=TEST_SAMPLE http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
<gml:boundedBy>
    <gml:Envelope srsName="EPSG:26912">
      <gml:lowerCorner>
        NaN NaN
      </gml:lowerCorner>
      <gml:upperCorner>
        NaN NaN
      </gml:upperCorner>
    </gml:Envelope>
  </gml:boundedBy>
</wfs:FeatureCollection>


texts: What means the gdal_translate message: "Too many command options"?

texts:
I have a little question, regarding the gdal_translate-utility:

I want to extract raster-files from database-stored gridfiles using the gdal_translate:

  gdal_translate -ot Float32 -of SAGA PG: "host='localhost' port='5432' dbname='$dbname' user='$dbowner' schema='grids' table='testrun' mode='2'" /output/...

But instead of delivering output it shows me "Too many command options." Then it lists usage and supported output formats. What does this message mean? I was searching the net already but found nothing which could be helpful.

Thanks for an answer in advance,
knutella


texts: I am using near feature tool in ArcGIS 10.1. However, I need to add another condition to it (not only the nearest). The feature (point) that i am looking for, is the nearest (point) and also doesnot cross another feature (line or area feature). For example, if the shortest distance line crosses another feature (line or area feature), then it has to start looking for another near feature (point) to which the shortest distance doesnot cross the feature (line or area). Is it possible using python?


texts:
I'm looking to accomplish a similar task to this thread: Splitting line at point positions using QGIS?

Using QGIS and the Processing and GRASS toolboxes available, I'm looking to split a vector layer representing a hydrolic network at points which represent junctions. I've been told v.distance exists in Processing as well, but I can't access it. Using GRASS v.distance, I get the error notice that not enough columns are available.


texts: Label Expressions to replace null values with text not in feature class

texts:
I need to label some hydrants (points) - some of these belong to the fire brigade and have a number assigned to them, some do not and thus do not have a number. I would like to write an expression along the lines of:

'Label Hydrant with number if number exists, if hydrant has null value label with 'WO'.'

I just can't figure it out.

The simplest way for me to do this would be to change the null values to WO but the data is not mine, and I do not want to make changes to it. I suppose I could make a copy and put the values in the copy, but ideally I would like an expression.

Thanks!


texts: How can we work with QGIS Server on Windows?

texts: I want to work with QGIS Server on Windows. I have not information about properties in this area. Where should I start?


texts: Changing value of point according to points' value inside its buffer using QGIS?

texts:
I have a layer with points and every point has a value.  I would like to create one buffer with a specific radius for each point.

In addition, I would like every point to acquire a new value the average of the points' values for the points which are inside the buffer.  If the point's buffer doesn't include any other point, the point's value will not change.

My professor told me about Spatial Moving Average and Floating Catchment Area but I have no idea how to use them.

I use mainly QGIS.


texts:
I am wondering what I would need to do in order to import a CSV file and then combined that with an existing layer.

I would like to know which boxes I need to check in order to make this combine with another layer.

If you could also walk me through the step by step procedures if we need to write any expressions in the field calculator that would be great.

If it helps I downloaded the original Layer with all the Data points on it and then Added a field and would like to re-add it into QGIS.


texts: Are there examples of high-level interfaces to NetCDF

texts:
The R raster package and GDAL both have independent wrappers around NetCDF-4* that treat an array variable as a "complete" dataset, i.e. they provide objects or interfaces that seamlessly give access to a variable inclusive with:

dimension variables
coordinate system metadata
value units
general other metadata

* I'm also interested in HDF4 and HDF5 extensions, but it seems less likely than for NetCDF-4.

I often see code that uses these at the lowest level, and I'm hoping to find more examples like those in raster and GDAL that have abstracted this away:

## psuedo code
load NetCDF
nc = open("file.nc")

v = getvar(nc, "somevariable")

x = getvar(nc, "somevariable_lon")
y = getvar(nc, "somevaraible_lat")

etc.

In R's raster and in GDAL it is much more like an extension of a traditional GIS raster, i.e. all the spatial context and metadata is present, and it extends to 3D and above. (There are problems with different conventions for rectilinear and curvilinear coordinates, but that's not my concern here).

Are there other open-source examples in widely used languages that provide the high-level interfaces like R's raster and GDAL?

Is there capacity in the NetCDF library itself for this higher level access? (I know how to use the lower level to build this, I want existing examples).

I'm also excluding uses of GDAL in Python, Perl, C# etc. - that is all just GDAL as far as I'm concerned. I'm also not interested in GIS interfaces that provide these via GDAL.

raster:
http://cran.r-project.org/web/packages/raster/index.html

GDAL: http://www.gdal.org

edit (2017):

Python's xarray is a good example


texts: How to merge polygons to get coordinates of final shape?

texts:
I have list of plots numbers that should be merged into multipolygon. As a result I have to get coordinates of final shape (as multipolygon).
That is may code:

    # work only on layer "plots"
name = "plots"
layer = QgsProject.instance().mapLayersByName(name)[0]
# fetch list of polygons
polyList = self.ptePlotList.toPlainText().replace(", ",",").split(',')
selected_fid = []
for poly in polyList:
      expr = QgsExpression("\"NUMBER\" = '" + poly + "'")
      it = layer.getFeatures(QgsFeatureRequest(expr))
      ids = [i.id() for i in it]
      selected_fid.append(ids)
selected_fid = [item for sublist in selected_fid for item in sublist]
# highlight selected polygons
layer.select(selected_fid)
# merge selected polygons and get coordinates, do not need to display
selection = layer.selectedFeatures()
geoms = QgsGeometry()
for feature in selection:
      geoms = geoms.combine(feature.geometry())
mergedPolys = geoms.buffer(0,5).asWkt()
print (mergedPolys)

but mergedPolys is empty. What I'm doing wrong ?


texts: ArcGIS 10.1 SP 1 install on Windows Server 2012 questions

texts:
I will be setting up ArcGIS server 10.1 SP 1soon and I have a few questions before I start this.  I have looked at the Web Adapter and the ArcGIS install guides but wanted to get some input.

1.Has anyone been installing on Windows Server 2012?  I will be installing ArcGIS 10.1 SP1 and the Web Adapter on a separate 2012 servers.  Have there been any problems with this that people have ran into?  Should I just stick with 2008 R2 or just go with 2012?

2.I see the ArcGIS server needs IE 7 min, and looks like it can go up to 9.  Will this software work with IE 10 though?

3.What ports am I going to need open between the Web Adapter server and the ArcGIS server, 6443?  Is 6080 just for the server to World if I am doing the web and GIS all on 1 server?

4.If I am using SSL and HTTPS over port 6443, should I be using certs?


texts: Finding lat/lon of corners of geotiff file

texts: I just recently managed to turn a GeoTIFF file into a .bmp file using gdal_translate.  Now I need to find the latitude/longitude coordinates of the corners.  I would prefer to use more of the gdal code or utilities, but that is not essential.  Does anyone have a good way to do this?


texts: NetCDF to Raster in ArcMap 10.3 - selecting Band Dimension stalls out

texts:
I am using the "Make NetCDF Raster layer" tool in ArcMap 10.3 and have had great success with relatively large files.  I usually select the Band Dimension (time) using the tool or sometimes after the raster is created under the layer properties NetCDF tab.  Either way it's been pretty fast and easy.

I'm now working with a much larger NetCDF file (2.3 GB; dimensions are 290 x 500 (x-y) and 1944 bands (time in hours)).  The NetCDF to Raster tool starts executing, and even says it's succeeded in the Results, taking about 22 minutes.  However, the output raster layer is empty, and the hourglass icon doesn't change (indicating completed).  It's now been > 6 hours and still nothing.

Is this file just too big for ArcMap to handle?  I was hoping to use the NetCDF_time_slice_export script tool to split my NetCDF up after getting it in as a raster, but I'm not sure that is possible.

Any suggestions are welcome.


texts: Table to point feature class using ModelBuilder of ArcMap

texts:
I am trying to convert a table in a fGDB to a point feature class using ArcMap 10.8. I need to do this in the ModelBuilder .
The table has about 8000 records, containing X and Y coordinates and many other attributes. I need to create points from the X and Y coordinates that also contain all the other attributes. I know the following procedure: "file -> add date -> add xy Date". But this does not provide a result that I can include in my model.
Is there any way I can create a point feature from xy data in ModelBuilder?