texts: Railways usage - How to assign values to part of existing lines?

texts:
I am conducting a project about rail accessibility/capacity, working with railways/stations downloaded in R with Overpass.

Manually I computed the number of trains for each lines and I would like to assign this values to specific segment of the SpatialLinesDataFrame.

For Instance, Zagreb-Osijek 20 trains, I want the line passing through this two stations to have 20 as value.

Obviously I don't want to create a new line, but to use the existing railway.

The problem is that the lines representing the railways have different starting and ending points from the station nodes.

Any suggestions to perform this process?

I thought that merging all the lines in one feature could help.

railways_merged <- gLineMerge(railways)

In this way I have only one feature, and I guess that I need to re-split them according to points (stations).


texts:
I use these options to create mbtiles for mapbox --accumulate-attribute to accumulate population from polygons, so I can get population with mapbox.queryRenderedFeatures in bbox
tippecanoe -o out_sum.mbtiles -P --coalesce-smallest-as-needed --detect-shared-borders --extend-zooms-if-still-dropping -z12 --accumulate-attribute=population_c:sum joined_output_z_F.geojson

but at small zoom levels like 1-2-3-4-5, some of the geometry is dropped.
Is it possible somehow to not drop geometry that accumulate-attribute?
I need the geometry(that accumulate-attribute population_c) visible at all zoom levels.


texts: What is the correct variable type for a single value in raster calculator, modelbuilder?

texts:
The equation includes several raster and several single values, and I want to do it in raster calculator, modelbuilder for the purpose of multiple calculations. In raster part, it is ok to select raster band, but how about the single values? I have selected any value under the function of create variable, but it does not work.

The expression is Cos("%Solar Zenith Angle%") * Cos("%Slope%") + Sin("%Solar Zenith Angle%") * Sin("%Slope%") * Cos("%Solar Azimuth Angle%" - "%Aspect%")

P.S. Slope and aspect are both raster, and Solar Zenith Angle and Solar Azimuth Angle are only single values.

It works, when I use raster calculator directly, so I think the problem probably comes from the variable type.

Any hint?
Thanks!


texts: Using World Files to Reference Raster Datasets in ArcMap

texts:
I am trying to georeference a simple raster image using a world file I was supplied with, however I am struggling to associate the world file with the image. I have read up on it, it seems that in the layer properties I should be able to 'edit' the spatial reference and I cannot. I have also specified in ArcMap options 'Use world file to define the coordinates of raster'. I am using ArcMap 10.4, I don't know if this will have made a difference.

Steps:
1) Checked the file names are the same (aside from the extension)
2) Customise - ArcMap Options - Raste Dataset - Check 'use word file to define  the coordinates of the raster'
3) Load Raster

What next?

Here is the contents of the PRJ and JGW files:


texts: SQL Query does not work in ModelBuilder (but works elsewhere)?

texts:
I am having issues getting a simple SQL query to work in ModelBuilder. It works fine when I use it with the standard "Select features by attributes tool", but ModelBuilder doesn't find the required fields.

"Sheet1$.Number" is not NULL AND Pipe.CLASS = 'IS17' AND "Sheet1$.CLASS" not like '17'

  ERROR 000230: Failed selecting with "Sheet1$.Number" is not NULL AND Pipe.CLASS = 'IS17' AND "Sheet1$.CLASS" not like '17'
  An expected Field was not found or could not be retrieved properly.
  An expected Field was not found or could not be retrieved properly. [Pipe]
  Failed to execute (Make Feature Layer).


texts: Spatial index for a raster column in a table in PostGIS

texts:
I have a raster table in PostGIS. When I run ST_SetValue(rast, geom , value) on it to change the pixel value at a point, it works very slowly.
Do you know how can I create spatial index for a raster column?
(I could not create it using PgAdmin interface with gist method)


texts: Reading PBF file or finding polygons of postcodes?

texts:
The issue: I want to find the polygons for the greek post codes.
I am not sure if this would help but I have just downloaded a PBF file from geofabrik.de.
I installed osmconvert and OsmAndMapCreator, initialy so that I can see what this file contains and then find a way to use, if in fact has something usefull there, since I am not even sure if it has the information I need.
Using OsmAndMapCreator I haven't found a way to load my downloaded file. I think it is connected to openStreetMap site.
Using osmconvert on the other hand I have managed to read the file, but I only get nodes. no polygons.
for example when I transform the file to csv

osmconvert greece-latest.osm.pbf --all-to-nodes  --csv-headline --max-objects=35000000  >test.csv

and try to find the boundaries (which i suppose I should get the country,city etch boundaries)
cat test.csv |grep boundary
all I get is

  node  1000000046421309    Territorial water boundary between Greece and Turkey
  node  1000000046421312    Territorial water boundary between Greece and Turkey
  node  1000000082288606    Territorial water boundary between Greece and Turkey

So, is there a way to see if postcode polygons are in the file?
If there are not, is there a way to get the polygons from openstreetMap or any other place?


texts: Alternatives to Offline Editing Plugin of QGIS?

texts:
I am looking for a way to edit my postgis layers offline (in field) using a Microsft Surface Tablet.

I have been testing the "Offline Editing Plugin" unsuccessfully. These are the issues I have experienced:

It does not work with WFS layers. When I try to convert the layer to offline QGIS crashes.
It does not work with MULTIGEOMETRIES
I get some warnings saying that it does not understand QDate fields (normal date fields in postgis).

Despite all these issues, I tried to use the plugin with a layer loaded directly from postgis and with a single geometries (No MULTIGEOMETRY). It turns the database offline Ok, I edit some features, but when I click on Sync data, the data is not updated in the database. Nothing happens in the database.

Due to this, I have discarded to use this plugin. Do you know any other software, architecture or whatever that could have a similar functionality, but working?


texts:
After creating a draw interaction:

var draw = new ol.interaction.Draw({
        source: vector.getSource(),
        type: /** @type {ol.geom.GeometryType} */ ('Polygon'),
      });

Is there a way to get its Style? In this case it would be the default Style.


texts: How much error should occur when using WGS84 to ECEF coordinates?

texts:
The method is not recommended for distances larger than 4km.  I am getting reasonable accuracy in terms of the east and north bearings, but for the up direction I am getting error upwards of 1km for distances of 110km (or a 1 degree difference in longitude).  The error seems to occur in the WGS84 to ECEF conversion, as the 1km error first appears there.  However, the result of the WGS84 to ECER conversion has been checked against numerous well known and presumably correct oracles such as the results of software like GeographicLib and the Sysense online conversion program.

For example using the WGS84 coordinates of (0,0,0) as a reference point I get the ECEF x,y,z coordinates of 6378137.00, 0.00, 0.00, this is correct.  So now using this as a reference point, I then use WGS84 coordinates (0,1,0) and get ECEF x,y,z coordinates of  (6377165.58, 111313.84, 0.00).  The difference between the two is (−971.42,111313.84,0.00).  This means that moving 1 degree east in the WGS84 coordinate system will give a -971.42 meter drop in x direction of the ECEF coordinate system (which is straight down relative to the initial position 0,0,0).  This cannot be possible because the curvature of the earth is about 8 inches per mile, and there are 68 miles in 110 km.   So the drop should be 544 inches or about -13.8 meters, not -971 meters.

Either the two oracles that the result has been checked against are incorrect, or the error is an inherent problem with the WGS84 conversion formula.


texts: Combining multiple custom QGIS plugins into one plugin?

texts:
I have built multiple individual QGIS plugins for my team using the QGIS Python plugin builder. Actually, my plugins are hosted on a private repository I made and the users can install them individually by adding the XML the repository to their Plugin manager settings.

Is there a simple way to make an entry in the repository XML or something else so that all my plugins are getting installed together?


texts: Does changing a .shp on ArcGIS desktop update the published ArcGIS Online map?

texts: I have a shapefile on ArcGIS for desktop that contains a few thousand points. I also recently published this shapefile to a new map on ArcGIS online. After working with the map for a while, I realized that some of the points are currently gone now, and there are some points that don't exist that I need to add. Can I somehow edit the shapefile on the desktop and then update the published map online? Or is the only way to change them both to individually edit each of them?


texts:
I own a (very) small appliance service company. It's just me. I cover roughly a 20 mile radius from a my office.
Customers can sign up on my website for service calls anywhere within my service area. I use Google's Distance Matrix API to reject any calls that are more than 20 miles or 30 minutes from my office.
Generally this works pretty well, however sometimes I'll get a call at the far east border of my area, followed by a call at the far west border, followed by another call at the far east border, which is much more wasted time and miles than I want to allow.
The question is "how do I detect and define 'the same area'" so I can refuse or reschedule the call to eliminate driving all the way across the territory more than once a day?
I haven't been able to figure out how to tell if I'm covering the same general area (for example, the northwest corner of the county) more than once a day.
Does anybody have any suggestions what I can do/look for and/or how I can define "same area" in a way that would be useful for routing?
Just to make things more complicated, the area is loaded with lakes and rivers, so a location that's 1000' away on the other side of some water geography might be an hour or more drive.
UPDATE - Where practice and theory collided
Initially, I tried a floating window over the "last location visited" but a small window excluded too many customers and a larger window let the area covered "wander" across the state if customers happened to show up in order of their coordinates.
5 or 6 customers that happened to have a similar latitude and sequentially increasing or decreasing longitude (or the other way around) could leave me driving across large portions of the state, which was exactly what I was trying to eliminate.
So I implemented @nr-aus idea and just split the state into quadrants, within a maximum radius. This means the area for any particular day is known and acceptably small and that the location test for new customers becomes a couple of latitude/longitude checks.
It's quite practical and lightweight and has been working well.


texts: Label inside top left corner of polygons using MapInfo Professional?

texts:
Does anyone know of a tool or procedure to automatically place a label in the corner of a polygon, rather than at the centroid?

I am using Mapinfo and know that I can manually place the label there, but I have a lot of polygons and would prefer to do it automatically, especially as the polygons change every so often.

Even better would be a way to do this in ArcGIS Desktop.


texts:
https://code.earthengine.google.com/51d66ca06d0b08dbb5a7f19c10817252

Using theese Lines wont result in a Image with the Line Painted on it

var Line =ee.FeatureCollection(geometry)
var Pimage = ee.Image(0).paint(Line,ee.Number(100),ee.Number(10000))
Map.addLayer(Pimage)

in the Documentation it says following:

ee.Image.paint(featureCollection, color, width)

this:image (Image):

The image on which the collection is painted.

featureCollection (FeatureCollection):

The collection painted onto the image.

color (Object, default: 0):

Either the name of a color property or a number.

width (Object, default: null):

Either the name of a line-width property or a number.

I do not understand why no line is painted in Pimage, since I have converted the ee.Geometry to a ee.FeatureCollection


texts: Fill sinks and Strahler Order error in QGIS

texts: I am working on a DEM of a watershed with a geotiff file size of 35GB. I am trying to execute watershed analysis tools and am experiencing errors in doing so. I attached a screen cast for Strahler Order error. The execution was finished but there was no output from it. Is it because of the big file size? I would be thankful if someone could shed a light on how I could resolve the issues.


texts:
Please forgive a very, very elementary question.

I am using a PostgreSQL 9.3.5 database, and one of the tables stores longitude and latitude.  (These are two separate columns.)  I would like to do efficient searches of the form "find all rows with (latitude,longitude) within 20 km of a given point".

I understand that this version of PostgreSQL has an "sp-gist" indexing method for addressing exactly this sort of situation.  However, I can't seem to find an example or tutorial of how I go about actually building the index.  In particular, are (latitude,longitude) distance handled natively, or am I somehow supposed to provide my own C code to implement the haversine formula?  Do I need to convert the (latitude,longitude) into some sort of joint spatial variable?


texts: PyQGIS: Zonal statistics not being written due to too many digits

texts:
I am using Python in QGIS to calculate zonal (sum) statistics for an array of raster layers using the following code:

polygonLayer  = QgsVectorLayer('E:/folder/' + 'gadm0.shp','','ogr')
rasterFilePath = "E:/folder" + r'/merged' + str(x) + '_R.tif'

polygonLayer.startEditing()
zoneStat = QgsZonalStatistics (polygonLayer, rasterFilePath, str(x), 1, QgsZonalStatistics.Sum)
zoneStat.calculateStatistics(None)
polygonLayer.commitChanges()

However, during the code I get some warnings like "value x not successfully written for feature z, possibly due to too larger number with respect to field width". I guess this is owing to the fact that shapefiles have a limit in the number of digits each field can contain.

When I open the attribute table of the output shapefile, I notice that:

1) a number of features have 0 as the output value, which does not make sense and thus I guess are those affected by the error.
2) successfully written features have way too many decimal digits than those I need.

How can I round the sum statistic to avoid data loss and useless decimal digits?


texts: Evaluate string in field calculator QGIS

texts: Is there any way to evaluate a string as an expression in QGIS.  For example, I want to evaluate the string in Column 'hukg' as an expression.


texts: Raster statistics for polygon areas PostGIS

texts:
I'm trying to get zonal raster statistics for each polygon area.
I have used the following method from this site:

SELECT bufid, (ST_SummaryStatsAgg(ST_Clip(rast, geom, true))).*
FROM rasttable, buftable
WHERE ST_Intersects(rast, geom)
GROUP BY bufid

This works great for a test area with 11 polygons.
However, when I try to run it on the same raster with a larger table of polygons (count of 190,000) I keep getting the following error:

  ERROR:  rt_raster_from_two_rasters: The two rasters provided do not have the same alignment
  CONTEXT:  PL/pgSQL function st_clip(raster,integer[],geometry,double precision[],boolean) line 8 at RETURN

How might I go about resolving this?

I'm using PostGIS 2.2.0 on PostgreSQL 9.4.4, Windows 7 64-bit.


texts: Restricting values to domain codes (beyond attribute table)?

texts:
Geodatabase fields can be assigned domains. ArcGIS Desktop uses the domains as lookup tables (a.k.a. validation tables).

When users edit data in the attribute table or attribute editing window, the values are restricted to the list of codes in the domains. This is a useful tool for enforcing data integrity at the application level.

However, it is easy for users to accidentally circumvent the domains and enter erroneous values by using the field calculator, ArcCatalog LOAD tool, geoprocessing tools, etc.. This is the main cause of errors in geodatabase tables at my organization.

Is there a way to use domains to restrict values that can be entered (beyond the attribute table) to truly force data integrity?


texts: Why would I need to use GeoServer's WFS-T capability while I can directly alter the PostGIS data directly?

texts:
I am trying to understand the use case of WFS-T to alter data coming from a database such as PostGIS. It is normally recommended to use the WFS-T capabilities of a Web Feature server only if the data is coming from a database store (Concurrent editing of data coming from shapefile published by GeoServer would be problematic for example).

Moreover, working from a client such as OpenLayers,  I can collect user input into GeoJSON for instance, handover the data to PHP, and save it directly into PostGIS.

My question therefore would be, when/where and why do I need to use the WFS-T capability of GeoServer?


texts: User name and password error during ArcGIS Server connection from ArcGIS Desktop

texts:
I get this error message, when I try to connect to our ArcGIS Server. I can reach that from browser with the same address. I get this error, when I try to connect from ArcMap, without any loading or any sign of real connecting.

  We were unable to connect to .... Error: Make sure you type in your user name and password correctly.

What can be the problem?

I'm not sure about my password or user name are wrong. According to my memories, an administrator user, who was picked up into the server's manager site with admin role should be able to log in. I tried with disabled firewalls too.


texts: Can you make QGIS decorations scale dependent?

texts: I have a QGIS project that includes a grid decoration. Is it possible to display this decoration only when zoomed in? I.e. can decorations in QGIS have scale dependent visibility?


texts: Display clipped to unknown boundary in ArcMap

texts:
I have a drawing in ArcMap. The display of all layers is clipped to an unknown boundary that appears to be made of overlapping circles; see image below; this is a screenshot of mapping visible in the drawing. In the image, the background map should be visible up the extents of the image (and thus a rectangle of mapping should be visible):

This applies to all data types; rasters, feature classes etc. I have tried changing all settings in "Clip Options" in "Data Frame Properties". I cannot find any possible settings to change in "Map Document Properties".
What are the possible causes of this clipped display?


texts:
I need to obtain information (road surface typology, max speed limit etc.) on all the roads that make my route, from A to B.
The inputs of my problem will be waypoints, obtained by GPS coordinates.
I need to extract road information based on these spatial way points.
I know it's possible to find these types of information from overpass turbo, but it has not helped me, since I don't know how to extract the route path given by the points.


texts: Reading existing layers from an existing project using PyQGIS

texts:
I am opening a QGIS project in a stand-alone Python script, but the existing mapLayers are not being found.  The script is below.  I am using PyQGIS 3.4 and I haven't been successful at using other solutions posted to the web for older versions of the API.

from qgis.core import *
from pprint import pprint

if __name__ == '__main__':
    import sys

    project_name = sys.argv[1]

    QgsApplication.setPrefixPath('/usr', True)
    QgsApplication.setThemeName('default')
    qgs = QgsApplication([], False)
    qgs.initQgis()
    print(qgs.showSettings())
    project = QgsProject.instance()
    project.read(project_name)
    pprint(project.mapLayers())
    qgs.exitQgis()

If I use the Python console in the QGIS application.  Then I can see this:

>>> from pprint import pprint
>>> QgsProject.instance().read('/home/janderso/code/data/audi/audi_debug/agg-545/audi_southbound.qgs')
True
>>> pprint(QgsProject.instance().mapLayers())
{'lane_marking_lateral_edges_041b777d_5e98_478d_96e6_ed3bf590a7a4': <qgis._core.QgsVectorLayer object at 0x7f710d7f4e58>,
'lane_marking_longitudinal_edges_be9d3f0b_ed83_4191_a0ee_1413209985b6': <qgis._core.QgsVectorLayer object at 0x7f710d7f4ee8>,
'lane_marking_nodes_29a65426_e9d8_4b8e_a74b_582ed3223e1d': <qgis._core.QgsVectorLayer object at 0x7f710d7f4f78>,
'road_boundary_lateral_edges_4f6c6b66_5d4c_4e3d_bb2c_d538215701ff': <qgis._core.QgsVectorLayer object at 0x7f7181f0e048>,
'road_boundary_longitudinal_edges_c0489433_be82_401a_88f4_8875e06025cc': <qgis._core.QgsVectorLayer object at 0x7f7181f0e0d8>,
'road_boundary_nodes_9fbd6dd6_69e8_420e_a126_e22cfd279df2': <qgis._core.QgsVectorLayer object at 0x7f7181f0e168>,
'sign_nodes_01febdb7_a31c_4f55_9f36_9ca71a077835': <qgis._core.QgsVectorLayer object at 0x7f7181f0e1f8>,
'sign_sign_edges_32c2531d_89bc_4046_8b2b_8d355920509a': <qgis._core.QgsVectorLayer object at 0x7f7181f0e288>}

Does anyone know the correct way to load an existing project using PyQGIS 3.4?


texts: Getting the projection of a LAS file using liblas?

texts:
I am using the python api of liblas to read/write las files. I would like to get the spatial reference information of the las file using liblas, but have been unable to figure it out.

Can anyone offer some advice on how to go about getting this information?


texts:
We have GeoNetwork installed in version 3.4.3. For some reason on the start page, there is no possibility to change the filtering of metadata based on INSPIRE themes.

We have INSPIRE Directive option enabled and INSPIRE thesaurus is added to the catalogue. What are we missing? How can we enable this functionality?


texts: Creating a BATCH with multiple INPUTS (variables) to run a workbench that intersects LINES AND POLYGONS

texts:
I'm trying to create a batch file to run my workbench.fmw in which I have two input files: 8 pipeline shapefiles and 96 communal areas shapefiles. The workbench allows cutting pipeline relative to the areas by using "Spatial Relator" Transformer.

I need to do batch processing in order to run this workbench several times with different inputs automatically (by double clicking on the .bat).

I tried to create my DOS script in a very complicated way so I come to you for help and suggestions and if anyone has ever created a similar batch?


texts: Loading shapefile to a specific table in postgis using ogr2ogr

texts:
I have already created a geometry table in postgis. I want to load several shapefiles into that table using ogr2ogr. But I didn't find the parameter which could specify that table.

ogr2ogr -append -f "PostgreSQL" PG:"dbname=db and so on" shapefile.shp

That's it. How can I add some parameter to specify the table?


texts: I am digitizing in ArcGIS Pro and suddenly the 'create structure and difference polygon' I had used earlier are grayed out and I can't access them again.


texts: MarkerCluster OR SuperCluster OR Vector OR

texts:
currently I'm working with many GeoJSON-Point-Objects. For my testing I use Leaflet with plugins (MarkerCluster, SuperCluster) and to compare them I tried OpenLayes with a cluster-vector solution.

I'm generating 30000 points with a nodeJS server and this delivers the points to the clients over a restAPI.

question: what is the best way to display more
than 30000 points in a browser without a mapserver?


texts: Offline entire planet in Java Desktop app

texts:
Sorry if this question has already been asked, but I need your help.

In program that we are developing we need to integrate maps. Here are the requirements that are put for integration:

Maps should be off-line, no access to online
Must be present all world, so the user can go anywhere as on website http://www.openstreetmap.org/
Maps must to be integrated(embedded) into the program(no third-party systems install)
Maps data should be in reasonable aisles
Maps need to work on the user's computer(not server configuration)

I'll add that we at this moment interested in the rapid integration of(a deadline) and again, our program is developed under the Java for desktop. Also we interested in free solutions.

As I understand it we need three solutions:

Database that can contain info about entire planet(from planet.osm file)
Lib that can render(on the fly) a piece of info from database into the tiles
Lib that shows renders results to user

Please help in choosing solutions/libraries and specify the way in which these things have already been made.


texts: Title label errors in infowindow in CartoDB - help?

texts: Sometimes, for no discernible reason, the title labels in my CartoDB infowindows start malfunctioning. They show the underscore in the title label, and I am not able to click on the "Change title labels" button (the A with the pencil). Any idea of what could be going on with this and how to fix it?


texts: Convert distance from meters to EPSG 3857

texts: I need to convert distance from meters to units of EPSG 3857 in PostGIS. This is important for me to calculate a radius for ST_DWITHIN as I'm dealing with geometries. How can I do this?


texts: Making chart with related table using ArcGIS Pro

texts:
I have a feature class with some locations of drills. I also have a table with data how deep the drill went and some more technical data in it.

I used relate on the feature class to connect the feature class with the table.
I am trying to make a linechart of the data of this table and I want to show this chart in the popup.

There is no option to select the fields of the table (related) as a X or Y axis.

Anyone know how to make a chart based on the related table?


texts:
I have saved group of layers which are reading from SDE DB as layer file. When I open the layer file in any other MXD it takes the same SDE data source.
How can I make the layer file without a data source?


texts: Generating lines from one polygon to group of polygons (one-to-many network) using QGIS

texts:
I have a polygon layer where every polygon has a unique "ID", and another polygon layer where those unique "ID"s can be found in a column in a one-to-many relation.
I would like to have QGIS draw lines from the polygons in the first layer to the polygons in the second layer in a way that reflects the one to many relation.
Any advice on tools?


texts: Column headers ignored when importing .xlsx file to postgres via QGIS DB Manager

texts:
I'm importing Excel (.xlsx) files to postgres using the QGIS DB Manager "Import Layer/File..." tool. I've already successfully uploaded very similar Excel files without issue, but for one file in particular, DB Manager is ignoring the column headers and instead populating them as the first row of data instead of as the headers.

I have already looked for extra spaces and special characters, and the headers are in good shape.


texts: I want to create a clean user interface for my students to use QGIS. I tried to remove the menu bar using the customization in setting but I am not able to remove everything, e.g. Database, plugin, etc. Is there a way to remove them?


texts: How to impose constraints on drop lists of two fields of the same layer?

texts:
Is it possible to impose with pyqgis3.0 some constraints on a drop list of a field attribute in order to govern the drop list attributes of another linked field?
To better clarify this issue, I attach the following schema:
When I select the Attribute 1.1 of the field 1, then the drop list of the field 2 should be limited to: Attribute 1.2, Attribute 2.2, Attribute 3.2.
Ect.


texts:
I am creating an app for tracking trains. When the train comes near to a railway station the system should identify it.

To do that I want to gather GPS coordinates near the station and save them.  Then when trains send GPS coordinates to the system, it will compare saved coordinates and identify the station name ...

How to find all longitude and latitudes (all gps coordinates)within the line from start point to end point?

If line is short (like 100 m long) i want to find all gps points in that line, so that I can save them in a database ..

thats why i want collect gps points in  that line


texts: How to measure the area of a hexagonal bin generated using ggplot2?

texts:
I've generated a number of plots using the geom_hex and stat_summary_hex functions in ggplot2, and I've adjusted the number of bins until the plot "looked right." That said, I'd like to know the area each bin represents to help interpret the plot.

Below is a simplified example of the set of example data and plot:

    library(dplyr)
    library(ggplot2)
    library(rgeos)
    library(maptools)
    library(rgdal)
    library(hexbin)

    crime <-read.csv("https://opendata.arcgis.com/datasets/6eaf3e9713de44d3aa103622d51053b5_9.csv", stringsAsFactors = F)

    homicides <- crime %>% dplyr::select(long = ï..X, lat = Y, OFFENSE) %>%
      filter(OFFENSE == "HOMICIDE")

     coordinates(homicides) <- ~long + lat
     homicides.df <- data.frame(homicides)

    homicides.hex <- ggplot() +
      geom_hex(data = homicides.df, aes(x=long, y =lat),
               fun = sum, bins = 20) +
      coord_equal()
    print(homicides.hex)

Is there a function or method to measure the area (in square miles or kilometers) of an individual hexagon when the x and y axes are latitude and longitude?


texts: Looking for a Java library for tracking positions, and querying regions

texts:
I have a deceptively simple requirement to create a software component to record and update some mobile position information, then be able to perform simple queries concerning mobiles in a particular area. When I say 'mobile' think aircraft and ocean-going ships. When I say 'simple queries' think 'Give me a list of all the mobiles in some region' or 'Give me the list of all mobiles within X km of mobile Y'. A simple, yet basically complete API for this component in Java might look like this:

import java.util.List;

public interface MobileTracker {

/** Tell us that this mobile should be tracked. */
public void logonMobile( int mobileId );

/** Tell us that this mobile should no longer be tracked. */
public void logoffMobile( int mobileId );

/** Update the current position of the logged on mobile */
public void updatePosition( int mobileId, double lat, double lng );

/** Finds all of the mobiles in the given region (circular or polygon). */
public List<Integer> getMobilesInRegion( /* Region region */ );

/** Finds all mobiles within radiusKm km of mobileId's current location */
public List<Integer> getMobilesNear( int mobileId, double radiusKm );

}

I'm looking for a Java library that will help me implement this. Preferably one that is not GPL, and doesn't have a million dependencies for a bunch of other functions that I don't need.

I don't work in the GIS domain, so as I looked here for help, I find myself in terminology-overload. It seems to me I should be looking for a spatial-database. This data doesn't actually need to be persistent here, so I would prefer something I could embed and run in-memory only. This question looked promising, but after a bit of digging on the HatBox site, I see it uses JTS which uses 2-D math, which I don't think is accurate enough for my needs, as I'm dealing with distances of 100s or 1000s of km (maybe this should be another question).

My question comes down to:

Am I right in assuming I need a spatial-database, or are there better solutions? Can you make some suggestions?


texts: Selecting an item with the mouse in MapBasic

texts:
This may seem relatively simple, but I'm pretty confused. I have written a mapbasic tool that does things to the table 'Selection' depending on user choices in a popup menu, but currently you need to use the select too to select an item and then run the program. I was wondering how I'd go about selecting an object with a click then it instantly running the main part of the tool?

Thanks,
Alan


texts: FeatureLayer Scale Settings in ArcGIS Silverlight API

texts:
Is there a way to translating between resolution and scale in the ArcGIS Silverlight API?

I am dynamically adding a Feature Layer to my map which should have a Min Scale value of 10,000.  The feature layer object has a setting for Minimum Resolution, but nothing for Scale values.

Does anyone know how I can set the Feature Layer to only show data below that scale?


texts: Calculate weighted centroid of a raster layer using GRASS?

texts:
I'm trying to calculate the weighted centroid of some raster layers I have, which are basically patches of positive values in a sea of null values.

I know I could convert thoses patches to binary data (1 inside a patch, null outside), vectorize them, calculate the centroid of each patch and then the "centroid of the centroids". But I really need to keep a coefficient corresponding to the value of each pixel in the original raster layer (I want that the core of each patch with higher values influence the location of the centroid).

It is not a problem in my case that the centroid will be outside existing patches, but I simply don't know how to perform this calculation from the raster layer using GRASS. Does any one know how to to this?


texts: To get 'NA' instead of '0' band value while exporting clipped image on Google Earth Engine

texts:
My problem is that I need to export images that are clipped using some shapefiles. But when I export the images, for the masked pixels I get 0 band value instead of NA. How can I get NA?
When I use Rodrigo E. Principe's batch download tool I get the exact format that I needed (na for masked pixels) but I could not use the tool since it is not made for a single image & I could not change the name of the export which is needed in my case.
var geometry = ee.Geometry.Polygon([
    [[23.42929490686074, 34.69481738849795],
     [46.36874803186074, 34.69481738849795],
     [46.36874803186074, 43.44349138830614],
     [23.42929490686074, 43.44349138830614]]
]);

// Import country boundaries feature collection.
var dataset = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017")

// Apply filter where country name equals turkey.
var turkeyBorder = dataset.filter(ee.Filter.eq('country_na', 'Turkey'));

var cloudyImage = ee.Image('COPERNICUS/S2_SR/20190427T083601_20190427T083603_T36SVJ');

var cloudBitMask = ee.Number(2).pow(10).int();
var cirrusBitMask = ee.Number(2).pow(11).int();

var qa = cloudyImage.select('QA60');

//Map.addLayer(qa.bitwiseAnd(cloudBitMask).neq(0), {}, 'clouds');
//Map.addLayer(qa.bitwiseAnd(cirrusBitMask).neq(0), {}, 'cirrus');

function maskS2clouds(image) {
    var qa = image.select('QA60');
    // Both flags should be set to zero, indicating clear conditions.
    var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
        qa.bitwiseAnd(cirrusBitMask).eq(0));
    return image.updateMask(mask);
}

var collection = ee.ImageCollection(ee.ImageCollection("COPERNICUS/S2")
    .filterDate("2019-05-01", "2019-05-06")
    .filterBounds(turkeyBorder))

Map.addLayer(collection, {bands: ['B4', 'B3', 'B2'], max: 2000}, 'collection');
Map.addLayer(table)
////Image bu
var min = collection.min().clip(table)
Map.addLayer(min, {bands: ['B4', 'B3', 'B2'], max: 2000}, 'min');
Map.centerObject(geometry, 5)
var bands = min.bandNames().getInfo();
for (var i = 0; i < bands.length; i++) {
    print(bands[i])
    Export.image.toDrive({
        image: min.select(bands[i]),
        description: 'may1-5_' + bands[i],
        folder: "may",
        scale: 100,
        region: geometry.getInfo(),
        maxPixels: 300000000000
    })
}

This is the link for my code https://code.earthengine.google.com/020cfd1f9e7a49ac325d0552a64f9f56


texts: Exported GeoPackage from R not drawing in ArcMap

texts:
I am trying to export multiple layers of point data from R as a single GeoPackage. The code works fine (included below), however when I try and load the resulting GeoPackage into ArcMap (10.6) I get the following pop-up:

Followed by this error:

I've not had this issue before when exporting GeoPackages from R and loading them in ArcMap and am not sure what is causing these problems (especially when it loads into QGIS fine). I have tried both writeOGR and st_write to export as a GeoPackage but encounter the same issue with both.

require("raster")
require("sf")

Data1 <- read.csv(sprintf("https://docs.google.com/uc?id=%s&export=download", "1RHAqtsPXgD5rH-ppIu7cLblHXyXzv_39"))
Data2 <- read.csv(sprintf("https://docs.google.com/uc?id=%s&export=download", "1Mzi6A6L_cXOrT9dTbm8gggGdtNQnZYkY"))
Data3 <- read.csv(sprintf("https://docs.google.com/uc?id=%s&export=download", "1ZzYiFIJxs2Coeb4kBNsWhmiF5NZDUhrL"))

Data1_Coords <- cbind(Easting = as.numeric(as.character(Data1$OSEAST1M)),Northing = as.numeric(as.character(Data1$OSNRTH1M)))
Data2_Coords <- cbind(Easting = as.numeric(as.character(Data2$OSEAST1M)),Northing = as.numeric(as.character(Data2$OSNRTH1M)))
Data3_Coords <- cbind(Easting = as.numeric(as.character(Data3$OSEAST1M)),Northing = as.numeric(as.character(Data3$OSNRTH1M)))

Data1_Points <- SpatialPointsDataFrame(Data1_Coords, data = Data1, proj4string = CRS("+init=epsg:27700"))
Data2_Points <- SpatialPointsDataFrame(Data2_Coords, data = Data2, proj4string = CRS("+init=epsg:27700"))
Data3_Points <- SpatialPointsDataFrame(Data3_Coords, data = Data3, proj4string = CRS("+init=epsg:27700"))

Data1_Points_ST <- st_as_sf(Data1_Points)
Data2_Points_ST <- st_as_sf(Data2_Points)
Data3_Points_ST <- st_as_sf(Data3_Points)

st_write(Data1_Points_ST, "DataGPKG.gpkg", "Data1")
st_write(Data2_Points_ST, "DataGPKG.gpkg", "Data2", update = TRUE)
st_write(Data3_Points_ST, "DataGPKG.gpkg", "Data3", update = TRUE)


texts: Selecting point with maximum value for each polygon

texts:
I have a point layer with two fields named "id" and "price", and a polygon layer with a field named "code_zone".

I need to get, for each polygon, the "id" of the point that has the highest value of "price", and this maximum value of the "price".

This is the beginning of my code:

for feature in layer_point.getFeatures():

   point_price = feature.attribute('price')
   point_geom = feature.geometry()

   for feat in layer_polygon.getFeatures():

      polygon_codeZone = feat.attribute('code_zone')
      polygon_maxValue = feat.attribute('max_value')
      polygon_geom = feat.geometry()

      if polygon_geom.contains(point_geom):
             ...

I tried with maximumValue(), but I get only the maximum value of the whole point layer and not the maximum value per polygon. What I am doing wrong?

As soon as I get these maximum values, I could fill them in my polygon layer with polygon_layer.changeAttributeValue() in the "max_value" field.


texts: readOnly attribute of ArcGIS FeatureTable

texts:
I am wondering the purpose and correct usage of the readOnly attribute of the FeatureTable widget in the ArcGIS Javascript API. The API defines it as a varible which "Indicates whether the data is editable via the widget." The API does not indicate it is a parameter that is set in the widget's constructor.

How is it used?

It would be great if there were a way to edit data within the FeatureTable object. However, I have tried both of the following implementations of code to no avail.

Nor have I found a code sample which uses this attribute.

Any ideas?

Option 1: Declare in the constructor despite API documentation

projectsTable = new FeatureTable({
    "featureLayer" : projects,
    "selectionMode": "single",
    "readOnly": false
}, 'grid');
projectsTable.startup();

Option 2: Manually declare after the FeatureTable is instantiated

projectsTable = new FeatureTable({
    "featureLayer" : projects,
    "selectionMode": "single"
}, 'grid');
projectsTable.readOnly = false;
projectsTable.startup();


texts:
I'm searching for the value from the field "demA_demA" which has the lowest elevation by having the value 8 from the field "outTable_2".

How can I do this in the following script.

import arcpy
import os
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
env.overwriteOutput = True

env.workspace = r"D:\Julia_T\projekt\aufgabe2\testdatenC"

inASCII = r"D:\Julia_T\projekt\aufgabe2\testdatenC\2007040_cloud_free.asc"
outRaster = inASCII + "_out.img"
rasterType = "INTEGER"
# Process: ASCII in Raster
arcpy.ASCIIToRaster_conversion(inASCII, outRaster, rasterType)
# Process: Projektion definieren
   arcpy.DefineProjection_management(outRaster,"GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHER OID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
# Process: Sample
outTable = r"D:\Julia_T\projekt\aufgabe2\ergebnisse3\outTable.dbf"
Sample(outRaster, outRaster, outTable, "NEAREST")
#works
# Join the table
# Table view
arcpy.MakeTableView_management(r"D:\Julia_T\projekt\aufgabe2\ergebnisse\outTable.dbf",     "out_tview") #works

joinField = "Y"
joinTable = r"D:\Julia_T\projekt\aufgabe2\ergebnisse3\demA.dbf"
arcpy.AddJoin_management("out_tview", joinField, joinTable, joinField, "KEEP_ALL")
# Persist the layer on disk
arcpy.CopyRows_management("out_tview",     r"D:\Julia_T\projekt\aufgabe2\ergebnisse3\joinTestC.dbf")
# Execute CreateFileGDB
arcpy.CreateFileGDB_management(r"D:\Julia_T\projekt\aufgabe2\result", "test1.gdb")     #works
outLocation = r"D:\Julia_T\projekt\aufgabe2\result\test.gdb"
# Execute TableToTable
arcpy.TableToTable_conversion(r"D:\Julia_T\projekt\aufgabe2\ergebnisse3\joinTestC.dbf",     outLocation, "jointest1") #works
# search for the lowest elevation with the value 8
with arcpy.da.SearchCursor (r"D:\Julia_T\projekt\aufgabe2\result\test1.gdb\jointest1",     ["outTable_2", "demA_demA"], """""demA_demA" = min("demA_demA")AND "outTable_2" = 8""") as     cursor:
    for row in cursor:
        print row [0], row [1]

Do I have to add the geodatabase because of an subquery? How should I write the subquery? I'm pretty new and working with ArcMap 10.2.1


texts: Debugging error 999999 from Python Script using feature layer?

texts:
I'm attempting run a calculate field for a feature layer that I'm creating in this tutorial. I can't seem to set it up to run though.                                                        ERROR 999999: Error executing function. Failed to execute (CalculateField)...     I want the expression to run and populate the new field in the feature layer,fc.

#!/usr/bin/python
# -*- coding: utf-8 -*-

# Import the mod
import arcpy
from arcpy import env

try:
    # Set up the environment
    env.workspace = \
        r"E:\GISTutorial_for_PyScripting\GISTPython\Data\FireDepartment.gdb"
    env.overwriteOutput = True

    # Create a list of all the box zone feature classes
    fcBoxZones = arcpy.ListFeatureClasses('FireBoxMap*')

    # fcBoxZones = arcpy.ListFeatureClasses("FireBoxMap_14")
    # Start a for statment to iterate through the files

    for fc in fcBoxZones:
        # Get the first file - it's stored in fc
        print fc
        break  # If you want the first otherwise it's the last

    # Add field to hold the results
    arcpy.AddField_management(fc, 'LaneMiles', 'LONG')

    # Select the street lanes
    path = r"E:\GISTutorial_for_PyScripting\GISTPython\Data\City of Oleander.gdb\Street_Centerlines"
    arcpy.MakeFeatureLayer_management(path, 'Lanes_lyr')

    # Replace a layer/table view name with a path to a dataset
    # (which can be a layer file) or create the layer/table view
    # within the script
    # The following inputs are layers or table views:
    # "FireBoxMap", "Street Centerlines"
    arcpy.SelectLayerByLocation_management(
        in_layer='Lanes_lyr',
        overlap_type='INTERSECT',
        select_features=fc,
        search_distance='',
        selection_type='NEW_SELECTION',
        invert_spatial_relationship='NOT_INVERT')
    laneCount = int(arcpy.GetCount_management('Lanes_lyr').getOutput(0))
    print str(laneCount)

    # Update field
    arcpy.CalculateField_management(in_table=fc, field='LaneMiles',
                                    expression=' [Lanes] * [SHAPE_Length]',
                                    expression_type='VB',
                                    code_block='')
    print 'updated LaneMiles field'

except Exception, e:
    print e


texts: Drawing point and path graphic shapes in print layout of QGIS?

texts:
At QGIS 3.4.4, I would like to freehand draw points and paths (i.e. zigzag lines) graphic shapes on my print layout.  However, I only see the option to draw rectangles, ellipses, and triangles, as shown in this screenshot:

How do I add points and paths?

Please note that I am not talking about digitizing these onto my map; rather, these will be simple graphics, drawn by hand over non-map portions of the layout, such as the legend, etc.


texts: Erase features from feature classes based on record ID?

texts:
This is the challenge:
I have two feature classes, which have the same number of polygon features, that when overlaid form concentric "rings." I need to erase the inner circle from the outer circle for each feature, matching them by their record ID. So, if both feature classes 1 and 2 have records A, B, and C, I want to erase the polygon of Record A in FC2 from Record A in FC1, RecordB in FC2 from Record B in FC1, etc.

ArcGIS' Erase function will only take the entirety of the feature class and erase from them - I can't find a way to link up records by their feature ID. I'm working with Arc10.1 with all licenses.

LATEST EDITS:
I've replaced the useless code from before with this revised snippet below, but now I need a bit of guidance - I can't find where in here that I am over-writing the output. It works! However the script leaves me with a single donut feature output. Can you help me spot my error?

#Test script: select record from a VTR and erase it from another, matching selected record

import os
import sys
import arcpy
from arcpy import env
import time

#Set overwrote output
arcpy.env.overwriteOutput = True
print arcpy.env.overwriteOutput

workspace = "C:\\BOEM\\pt_buff3\\"

fileList = []
buffersd4 = []
buffersd3 = []
recordID3 = []
recordID4 = []

buffs = workspace+"buffer.gdb"
directorypath = buffs

#create new results_feature class

for dirpath, dirnames, buffs in arcpy.da.Walk(directorypath,
                                              datatype="FeatureClass",
                                              type="Polygon"):
for buff in buffs:

    if buff.endswith(str("d4")):
       # buffersd4.append(buff)
        #print buffersd4
       print buff
       buffbits4=buff.split('_')
       print buffbits4[2]
       buffersd4.append(buff)
       field = ['APR16_ID']
       APR16_ID = [row[0] for row in arcpy.da.SearchCursor(workspace+"buffer.gdb\\"+buff,(field))]
       print APR16_ID
       for ID4 in APR16_ID:
           ID4split = ID4.split('_')
           selectID4 = ID4split[1]
           #print selectID4

    if buff.endswith(str("d3")):
       # buffersd4.append(buff)
        #print buffersd4
       print buff
       buffbits3=buff.split('_')
       print buffbits3[2]
       buffersd3.append(buff)
       field = ['APR16_ID']
       APR16_ID = [row[0] for row in arcpy.da.SearchCursor(workspace+"buffer.gdb\\"+buff,(field))]
       print APR16_ID
       for ID3 in APR16_ID:
           ID3split = ID3.split('_')
           selectID3 = ID3split[1]
          # print selectID3
print selectID4
print selectID3
if selectID4 == selectID3:
    print "yay they match!"

    #make feature layer of feature d4 for FCs in buffersd4:
    for fc in buffersd4:
        in_feat=workspace+"buffer.gdb\\"+fc
        out_layer4= fc+"_lyr"
        arcpy.MakeFeatureLayer_management(in_feat,out_layer4)
        print arcpy.GetMessages()
    #select record in fc_layer
        #where_clause4=""" "APR16_ID" = """ + str(selectID4)
        where_clause4=""" "APR16_ID" = """ + "'"+ID4+"'"
        arcpy.SelectLayerByAttribute_management(out_layer4,"NEW_SELECTION",where_clause4)
        print arcpy.GetMessages()

    #make feature layer of feature d3
    for fc in buffersd3:
        in_feat=workspace+"buffer.gdb\\"+fc
        out_layer3= fc+"_lyr"
        arcpy.MakeFeatureLayer_management(in_feat,out_layer3)
        print arcpy.GetMessages()
    #select record in fc_layer
        #where_clause3=""" "APR16_ID" = """ + str(selectID3)
        where_clause3=""" "APR16_ID" = """ + "'"+ID3+"'"
        arcpy.SelectLayerByAttribute_management(out_layer3,"NEW_SELECTION",where_clause3)
        print arcpy.GetMessages()

 #execute Erase function(erase d3 from d4)
        out_feat = "C:\\BOEM\\Donuts\\donuts_output.gdb\\"+"donut_d4_"+ID4
        arcpy.Erase_analysis(out_layer4, out_layer3,out_feat,"")
        print arcpy.GetMessages()

    print "done rolling through records in this feature class."

    print "done going through the buffers in the fGDB"

    #celebrate

            `


texts: In ArcGIS 10, how to use Append tool to "transfer" the fields of a shapefile to another?

texts:
I am at the end of a small project: conversion of a shapefile with a large number of features, therefore huge dbf table. The initial shapefile, after processing, ends up with fields of A, B, C, D, E, F etc. And now I want certain fields of them, about 20. But I want their name to be different, say a, b, c, d etc.

I know I can manually alter the field name but it's gonna be very time consuming and sort of silly.

My approach in mind is I will try to create a feature class of a shapefile with fields a, b, c, d but no features, therefore empty, and then use Append to kind of "assign" the A, B, C, D fields in the first shapefile to the a, b, c, d fields of the second shapefile instead of changing names of 30 fields one by one.

Is this doable? If so, how shall I create an empty shapefile with the fields that I need? Is there any better solution to this situation?


texts: ArcGIS 10.2 Query Layer on SQL Server performance

texts:
I'm using a Query Layer on SQL Server in ArcMap. The Query Layer executes instantly in in SQL Server but but takes so long to draw in ArcMap that the system appears unresponsive for about 10 minutes or longer. During the ArcMap draw the one of the CPU's is maxed out on the SQL Server process.

My Query is the STIntersects of a buffer on a line feature (Shannon) against a polygon feature class (Townlands) , as follows;

SELECT TOWNLANDS.TL_ID,TOWNLANDS.Shape FROM dbo.TOWNLANDS as townlands
with(index(FDO_Shape))
JOIN dbo.Shannon on townlands.Shape.STIntersects
(Shannon.Shape.STBuffer(2.0))=1

The query returns 186 rows instantly. These can be draw in the SQL Server Management Studio Spatial pane without a problem

When I build a Query Layer in ArcMap with exactly the same syntax the system becomes unresponsive but does draw eventually. It looks as though, perhaps, ArcMap s not using the spatial index or is doing so differenlt from SQL Server causing an inefficient query on SQL Server that takes an age to return.

Can anyone advise on a remedy?

Thanks

ArcGIS Desktop: 10.2
ArcSDE: 10.2
RDBMS: Database and version: SQL Server 2008
OS: Windows Server


texts: Programming language of QGIS field calculator?

texts:
For the functions of the field calculator in QGIS, in which programming language does it have the most parallels?

In ArcGIS field calculator you have to decide between VBScript, JavaScript and Python.


texts: Where is the Bathymetry toolbox in ArcMap

texts: I can't find the Bathymetry toolbox on my ArcMap 10.4. I went to the customize toolbar but there is no sign of it. I searched in but nothing shows up.


texts: Converting SWEREF99TM coordinates to the "standard" coordinate system (longitude and latitude)

texts:
I recently asked a question about a coordinate system in a file, and found out that it used SWEREF99TM, but I can't figure out how to convert those coordinates to the "geographic coordinate system". I've found a bit about it on lantmateriet.se. It says that:

  It is an advantage for many different applications at the national
  level if the whole Sweden can be mapped seamlessly, i.e. on one
  projection with position given in one co-ordinate system. A rigid
  application of, for example, UTM in Sweden would require a division
  into three zones.

  The definition of a dedicated projection based on a meridian 15° E of
  Greenwich, with a scale factor of 0.9996, with the Equator as the
  origin for N co-ordinates and with E co-ordinates measured from the
  central meridian and increasing eastwards (the central meridian is
  assigned a value of 500 000 m to avoid negative values) will give
  consistency with UTM in the major part of the country.

I still can't really figure out how to convert it, but I would be glad if someone could help me with it. I got a suggestion to use the java library proj4j (because I'm going to convert it using java, and then use the coordinates in an android app).
But I can't figure out how to use that either, and I also want to do it myself. I don't really like libraries. I want to do stuff myself. Don't know why. It just is like that :)

What does the scale factor mean, and what unit is the coordinates measured with, and so on.


texts: How to analyze Vegetation Index using Sentinel-1 SAR Imagery?

texts: I am using the ESA open scientific data hub for Sentinel 1 SAR imagery (Interferometric Wide swath SLC VV polarization) and would like to know how to perform vegetation health analysis from this data?  What tools are available for this using ESA's SNAP open source software?


texts: Determining hardware requirements for mosaic dataset operation to work properly?

texts:
I've been building mosaic datasets in ArcGIS, but have been having mixed results depending largely on the number of cores on the computer being used to build the mosaic dataset. For example, on a computer with 2 cores the operation fails, but with 4 cores it works.

I've found relatively little on Esri's website relating the capabilities of the mosaic dataset to the hardware of the computer from which it is being created. This leaves us with only trial and error to figure out the hardware requirements for a mosaic dataset operation to work properly.

I'm wondering if anyone has created a guide or knows of a guide that can be used to estimate the kinds of/size of operations that can be undertaken on a computer of a given number of cores, RAM, etc?

Specifically, it is the "add rasters to mosaic" operation that is failing.


texts: Compute the number of shortest paths that use a link

texts:
I am currently computing shortest path algorithms between two sets of points (origin set and destination set) and a road network, using the plugin QNEAT3. It works and gives me the distances between each point using the road network.
Still, my goal would be to calculate the number of shortest paths that go through a link, to calculate the traffic flow on the network.
I tried to use AequilibraE plugin but I cannot make it work with my layers for now. Do you have hints on what I could do?


texts:
I want to plot two countries:

library('maps')

map(database='world',regions='germany')
map(database='world',regions='burkina faso')

But "burkina farso" is clipped by the x-region of "germany".
How can I get a complete map of "burkina faso"?


texts: Raster Conversion set output path [GDAL]

texts:
I planned to use gdal_translate to convert raster from uncompressed RGB geotiff to compressed(JPEG compression) YCBCR geotiff as a batch process. So I have created a .bat file using this code

for %%i in (*.tif) do gdal_translate -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -co TILED=YES %%i %%.tif

The code works fine but the problem is output rasters are created in the same path of input rasters. For example If input raster name is RGB.tif, after processing .bat file output created in the same folder with the name of RGB.tif.tif.

Is there any possibility to store output tiff in a separate directory? I have kept my input rasters in directory D:\work\RGB and I want to store output rasters in D:\work\YCBCR with the same name of Input.


texts: QGIS crashes when starting Python, unless opened in administrator mode

texts:
I am using QGIS 3.14.16 and I believe, I have somehow corrupted my QGIS/Python interaction.
My Problem:
When I open QGIS and the status information reads "Starting Python", the application just closes. No error message, it is just gone. When I start in administrator mode however, it all works fine.
How I ruined the setup
I had it running just fine after the initial installation. I then changed the PYTHONPATH and PYTHONHOME variables in QGIS to try to work on a virtual environment on my system. After restarting QGIS, it would not start again. I have not added any System or User Variables in the Windows settings.
How I tried to fix it
I have uninstalled and deleted all files associated with QGIS twice now, but whenever I reinstall the problem is the same. QGIS is now installed in C:\Program Files\QGIS.
I checked the file \QGIS\bin\qgis-bin.env and PYTHONHOME and PYTHONPATH look fine to me:
PYTHONHOME=C:\PROGRA~1\QGIS\apps\Python37
PYTHONPATH=C:\PROGRA~1\QGIS\apps\Python37;C:\PROGRA~1\QGIS\apps\Python37\Scripts

The content in \QGIS\bin\p3_env.bat also looks fine to me:
SET PYTHONHOME=%OSGEO4W_ROOT%\apps\Python37
SET PYTHONPATH=%PYTHONHOME%;%PYTHONHOME%\Scripts
PATH %PYTHONPATH%;%PATH%

My research
I found this unresolved Bug, which describes someone messing with PYTHONPATH and PYTHONHOME as well and consequently resulting in QGIS crashing. I believe he changed the variables not within QGIS (what I did), but within his User Environmental Variables. After deleting them, it worked for him again.
I found a similar bug here where the solution again was to delete environment variables. Again, I did not set any windows environment variables.
As I said, starting in administrator mode is possible, but this is not my preferred option. I want to make sure that my setup is well defined. Additionally I have encountered problems when working with "Actions", which might be related to this problem.
Does anyone have ideas, what might cause the troubles? Any suggestions on what to try?


texts: QGIS KDE (Heatmap) TIF file --> Open in ArcGIS in order to use Kernel Interpolation With Barriers

texts:
So I've created KDE .tif files in QGIS, but there are currently no plugins available in QGIS, which takes into account barriers when creating KDE density rasters for species.

I can bring the .tif file into ArcGIS 10.4 and it is visible, however when I open the Kernel Interpolation With Barriers Tool in ArcGIS, the .tif file is not visible in the drop down menu [INPUT FEATURES].

I am wondering can I convert the QGIS .tif file into a file that is usable in ArcGIS? and can then be used with the Kernel Interpolation With Barriers Tool?


texts: Serialize python list containing geometry

texts:
In simple case if list doesn't contains geometry all works fine:

result = [dict(f.fields, id=f.id, label=f.label) for f in features]
json.dumps(result)

But if result list contains geometry (Shapely) I get error object at 0xb4a87bcc> is not JSON serializable:

import geojson
result = [dict(f.fields, id=f.id, label=f.label, geom=f.geom) for f in features]
geojson.dumps(result)

If I try to serialize not list but one dict - all works fine:

geojson.dumps(result[0])

Is it possible to serialize list of dicts?


texts: Applying a single leaflet search control to search multiple layers?

texts:
I am attempting to add a single search button to look through geoJSON layers. They will look through the field "name". Unfortunately, my method of creating a layer group is not working. What am I missing? Or how do I group together geoJSON files to make this work?

var searchControl = new L.Control.Search({
      layer: [syriaLayer, countryBoundary],
      zoom: 7.5,
      propertyName: 'Sheet_Numb',
      circleLocation: false
  });


texts: Google Maps - Use fitBounds on KML polygon fusion table layer

texts:
I have seen example here how to use map.fitBounds on fusiontablelayer
here and here , but I think it works with markers table only. Its impossible to parse the lat and long when the vertex coordinates are in the geometry column as in a kml polygon layer.

Is there a solution to use fitBounds or other function to zoom to the extent of a kml polygon fusion table layer?


texts: Earth Engine: Mapping over a feature collection of ROIs

texts:
I have a function that performs an operation over a specific ROI. Now, I want to apply the same function over a collection of ROIs by mapping the function over a feature collection.
Haven't found any feasible solution anywhere.

Here's what I am trying:

var myfunc = function(feature) {
  var geometry = feature.geometry();

  var L7 = ee.ImageCollection("LANDSAT/LE07/C01/T1_TOA")
      .filterDate('2009-10-01','2010-12-01')
      .filterBounds(geometry)

  var Convert_scale = function(img) {
    var I = ee.Image(img)
    var correct_scale = I.select(['B1','B2','B3','B4','B5','B7']).multiply(0.0001)
    return I.addBands(correct_scale,['B1','B2','B3','B4','B5','B7'],true)
  }

  L7 = L7.map(Convert_scale)

   return L7;
}

var fc = ee.FeatureCollection("users/climateClass/strips_bross_allin");
var L7coll = ee.ImageCollection(fc.map(myfunc))
print(L7coll) //.first())
Map.addLayer((L7coll).mosaic())

The resulting image collection has zero bands and adding first image to map gives a Layer error

Here's my code.


texts:
Recently QGIS (currently using 3.10 A Coruna) hangs/stalls on startup.  When I start the program, the splash screen appears and it runs through all the preliminary installation processes (javascript, plugins, etc.), then just hangs and doesn't actually open the QGIS program.  The splash screen says "QGIS Ready!", but then, the program doesn't open.
The only way I have managed to solve this is to uninstall the program, reboot, reinstall, and then start the program.  It starts correctly the next time, but then when I close it and shut down the computer, the next time I boot the computer and start QGIS, it just hangs like this.
I hate having to reinstall QGIS all the time.  Any suggestions about how to correct this problem without having to reinstall?
Thanks for your help.


texts: How to use creation_time and forecast_time in google earth engine to get forecast data?

texts:
I want to get the forecast rainfall for let say three days after this day of writing using the data from Google Earth Engine's NOAA GFS (Global Forecast System): 384-hour predicted atmosphere data. However, I am having a difficulty in utilizing the "creation_time" and "forecast_time" in the codes. I was told that this is different from the .filterDate command.

//Load input imagery.
var collection = ee.ImageCollection('NOAA/GFS0P25')
  .select("total_precipitation_surface")
  .filterDate("2018-02-15", "2018-02-19");

// Define a region of interest as a buffer around a point.
var geom = ee.Geometry.Point(124.265, 7.8675);

// Center map
Map.centerObject(geom, 14);

//Create time series and print the chart.
print(ui.Chart.image.series(collection, geom, ee.Reducer.mean(), 250));

//Create day of year (DOY) and print the chart.
print(ui.Chart.image.doySeries(collection, geom, ee.Reducer.mean(), 350));


texts: MapProxy zoom level based on the area of the MBR

texts:
I am working with MapProxy and I'm trying to determine the zoom level that I should use based on the area of the a rectangle (the MBR of a polygon).

How can I do something like that?


texts: Importing extension data from GPX files?

texts:
I have a Garmin GPS that does handy things like waypoint averaging which is very useful in my situation where much of my data collection happens under heavy canopy in deep gullies which leads to inaccurate fixes.  I often take several samples on different trips to get accurate positions.

The Garmin gpx files have a samples attribute in the extensions and I would like to be able to get this loaded into QGIS so I can have a script update the geometry if the sample number has changed.

here is a sample point with extension data

<wpt lat="-36.599535" lon="174.887440">
  <ele>2499999890505881600000000.000000</ele><time>2017-04-27T00:54:16Z</time><name>TJ B2 1</name>
  <cmt>Track junction in B2</cmt>
  <sym>Flag, Blue</sym>
  <type>Track</type>
  <extensions><wptx1:WaypointExtension><wptx1:Samples>4</wptx1:Samples></wptx1:WaypointExtension></extensions>
</wpt>

to load the file I am using the "add GPX layer" -- the old gps tool.

the GPS exposes the files as a USB file system so I am just reading them as files.  I can preprocess the files to, say put the sample data into the 'cmt' field but I would rather have a solution in QGIS because I want others to use this.


texts: Adding field values to a QComboBox in PyQGIS

texts:
I have a layer say 'Cities'. There is a column called 'NAME_1' in the layer. My GUI has a QComboBox in which i want the values of column 'NAME_1' alone to be displayed in the drop down of the combo box.

Can anyone help me in providing snippet in Python for getting a value from a particular column and displaying it in the combo box?


texts:
I am quite new to GIS so I am sorry for my "primitive" vocabulary.
In QGIS, I have a set of points falling in specific coordinates in Africa and Europe.
I would like to select only those points which are within the continents (thus not in the sea).

Probably the figure is more explanative :

I have tried to use the intersection feature but it is not working as the only result I am obtaining is the only point which is actually falling perfectly onto a line (see image. I have made the point bigger just to visualize it better).

Can you help me with this?

.png


texts: I have a raster file with pixel size 30 x 30. At the upper border there is half row with "no data" pixels. I want to use the mode function to replace them. I have only thought r.neighbors (from grass) with a moving window of 2x2 (just to take the values of the previous row). I haven't try this for filling areas on borders so my question, is it safe to use this function for borders where I have only one row as neighbors? (values are categorical of land cover)


texts: Over writing one large DEM over another in QGIS or GDAL

texts:
I want to over write/mosaic one large DEM (around 3gb file, Geotiff) over another. What is the best way to do this?

Both cover roughly the same region, one is based on 2m res LiDAR, the other 1m res LiDAR. Both data sets have large areas of no data (-9999), and both have a 1m cell size and the same projection. The resulting DEM should have the better quality LiDAR written over the lower quality where it exists.

I have tried using gdal_merge but this does not like the no data switch, takes too much memory (I do have 32gb of Ram). I have tried cutting down the files into smaller segments, but even at two 300mb files it still will not work.

I have tried using gdalwarp (-overwrite) but this simply overwrites the entire file giving me a copy of the 1m res data with no 2m beneath it despite using the no data switch.

I am open to different software, but the resulting file has to be converted into an ASC to be read into a model.


texts: How to SortBy a GetFeature WFS request by FeatureId (fid) in XML?

texts:
This is my XML for a GetFeature request. I want to Sort By FeatureId (fid) but this doesn't work. I can't find an example of how to format this XML correctly:

(SEE THE ogc:SortBy SECTION)

<?xml version="1.0" encoding="UTF-8"?>
<GetFeature xmlns="http://www.opengis.net/wfs" xmlns:gml="http://www.opengis.net/gml" xmlns:myws="http://mysrvr.com/solution/myws" xmlns:ogc="http://www.opengis.net/ogc" version="1.1.0" service="WFS" maxFeatures="100" outputFormat="JSON">
   <Query typeName="myws:mylayer">
      <ogc:SortBy>
         <ogc:FeatureId>fid</ogc:FeatureId>
         <ogc:SortOrder>ASC</ogc:SortOrder>
      </ogc:SortBy>
   </Query>
</GetFeature>


texts: It was a .h5 file. I converted the datasets to .tif format using gdal. Now I want to assign projection to those files but the prescribed projection of EPSG:9805 is not found in gdal. I tried other projections like Mercator (variant A) and web mercator, but they dont seem to fit when I overlay a shapefile of the same area on it. Can anyone tell if I am going wrong anywhere or suggest some other alternative.


texts:
I'm using the CartoDB static Maps API.  I need to change the font size of the street names.  E.g. in the image below the street name "Stanhope Road" needs to be bigger.

Here's what my Json config file looks like so far.  I've just been guessing how to set the cartocss property, and so far it hasn't worked.

{
  "version": "1.0.1",
  "name": "map",
  "layers": [
    {
        "type": "http",
        "options": {
            "urlTemplate": "http://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png",
            "cartocss_version": "2.3.0",
            "cartocss": "#road_label {  text-size: 24;  }"
        }
    },


texts:
We need to build a Web based GIS application for Raster data computation (analysis). All the maps are in form of raster images.

We are failing from a week to find a proper tool which allows fetching, viewing and computing of raster data easily. Let me mention here, We need to work on open source tools.

Can you please guide me to proper direction?

I have checked WCS on web. But frankly I am very unsure how to proceed further. A step wise guide would seriously help me a lot.


texts:
In the past I have used OpenLayers in order to make WFS requests to geoserver and show features on a map. I have also used WFS rules in order to filter the data based on specific values.
Now I am using Leaflet and I wonder whats the way to do something similar. Obviously Leaflet doesn't support many of the functionalities OpenLayers has but is there a way to show filtered features on the map?
Should I just bring the whole json through ajax jquery and then get the data I want from the jquery and show them?
Would it be an appropriate way?


texts: Problem with precision when convert from WKT to shapefile

texts:
I have a problem with converting from WKT to shapefile in pygdal. In wkt I have fields with precision (for example real 7.2). When I create a shapefile I have field numeric with no precision (numeric 7.0). I can't find problem.

Do you have any suggestions?


texts: Merging lines where point feature class exists in ArcGIS Desktop?

texts: I have a polyline feature class and a point feature class.  Two lines from the polyline feature class come together (in most cases they are snapped together) where a point exists in the point feature class.  Where two lines come together in this way, I need to merge those two lines.  Since the amount I of data I have is pretty tremendous, I'm wondering if there's a model or a python script that will help me to accomplish this much faster than doing it by hand.


texts:
My System (Ubuntu) previously had GDAL 1.11.2. We have an issue with using ogr2ogr converting KML to SHP in that the attributes were missing. In researching this, it seems that the LIBKML driver should fix this. I have tried researching various questions and tutorials on getting the LIBKML driver to work but I never seemed to get it working.

I have seen in other places that the LIBKML driver comes with GDAL in version 1.9.2 onwards. I got GDAL 2.1.3 on my system now, (verified with gdal-config --version), however when I run ogrinfo --formats, it is still only showing back KML and not both KML and LIBKML.

Is the LIBKML driver supposed to come with GDAL 2.1.3 out of the box?

Update:

As per request, here is the output from apt-cache policy gdal-bin

gdal-bin:
  Installed: 2.1.0+dfsg-1~trusty2
  Candidate: 2.1.0+dfsg-1~trusty2
  Version table:
*** 2.1.0+dfsg-1~trusty2 0
       100 /var/lib/dpkg/status
    1.10.1+dfsg-5ubuntu1 0
       500 http://mirror.rackspace.com/ubuntu/ trusty/universe amd64 Packages

Update with Solution
I downloaded gdal in a weird way, as pointed out by @AndreJ
Everytime I would uninstall everything, then reinstall the package, it still kept coming up with version 2.1.3 which isn't available through the PPA, which apparently doesn't have LIBKML.

What I had to do was manually locate gdal and locate kml and rm any and all files and instances found. After everything was gone, I could install gdal again and version 2.1.0 was installed with the LIBKML driver.


texts: Removing a polygon-defined area from a multiband raster: A landscape classification inquiry

texts:
I am trying to exclude multiple lakes and reservoirs from a Supervised landscape classification but at the moment I have not found a way to do this. So far I have used the following workflows:
1.Draw a polygon -> Clip -> Extract by mask) generates the opposite result.
2.Draw a polygon -> "Raster to Polygon" -> Reclassify -> Set Null, but the reclassification removes RGB symbology
Details:
I am using an Multiband High resolution Raster image
I am working on the latest version of ARCmap for Desktop
I would like to retain the RGB Composite symbology


texts: Symmetrical difference between two layers using PostGIS

texts:
This question is related to Seeking PostGIS equivalent to GeoPandas union overlay?

I have two tables in PostGIS, each with two features(rows) and a rectangular geometry (polygon).

I would like to obtain the Symmetrical Difference from those two layers so that the result would be:

the resulting table has the attributes from the first table (df1), or the second (df2)

However when I try to use ST_SymDifference in PostGIS I get a different result:

CREATE TABLE test.symdifference_v04 AS
SELECT
    test.df1.df1,
    test.df2.df2,
    ST_SymDifference(test.df1.geom,test.df2.geom)
FROM test.df1, test.df2

It seems that the ST_SymDifference is applied on a polygon to polygon basis i.e. df1 row 1 vs df2 row 1
df1 row 1 vs df2 row 2
df1 row 2 vs df2 row 1
df1 row 2 vs df2 row 2

with Multipolygons as a result.

How would I obtain the first mentioned result with the 6 polygons (fancy colors) instead of the 4 multiploygons in PostGIS?


texts:
am following this tuto to create a map on web with geoserver, GeoEXT, EXTJs, and openlayers      here's the linki downloaded ext 4.2.0, openlayers 2.12 and GeoEXT 1.1, when i try to show the map : //geo.ancfcc/geoserver/www/search/map.html , all i get is a blank page.
i don't know what's the error (i don't know where to find error log)
here's my map.html:

   <!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01//EN">
    <html>
     <head>
       <title>
        A map
       </title>
<script src="../ext/ext-all-debug.js" type="text/javascript"></script>
<script src="../ext/bootstrap.js" type="text/javascript"></script>
<script src="../ext/ext-all.js" type="text/javascript"></script>
<link rel="stylesheet" type="text/css" href="../ext/resources/css/ext-all.css">
<script src="../openlayers/OpenLayers.js" type="text/javascript"></script>
<script src="../GeoExt/lib/GeoExt.js" type="text/javascript"></script>
<link rel="stylesheet" type="text/css" href="../GeoExt/resources/css/geoext-all-debug.css">

<script type="text/javascript">
     Ext.onReady(function() {
      var map = new OpenLayers.Map();
      var layer = new OpenLayers.Layer.WMS(
        "Blue Marble",
        "http://geo.ancfcc/geoserver/wms",
        {layers: "ancfcc:titres_salemed_valid"}
      );
      map.addLayer(layer);

      new GeoExt.MapPanel({
          renderTo: 'gxmap',
          height: 400,
          width: 600,
          map: map,
          title: 'A Simple GeoExt Map'
      });
  });
</script>


texts: D3 Map + TopoJSON (mapshaper) showing empty

texts:
Just started learn D3.js and TopoJSON. I had shapefile files and converted using mapshaper to TopoJSON, but when I tried on my local server my map showing empty page without any errors. I don't have direction what's going on.

Here my HTML got from a tutorial:

<!DOCTYPE html>
<meta charset="utf-8">

<!-- Set a style for our data -->
  <style>
  path {
    stroke: red;
    stroke-width: 0.5px;
    fill: grey;
  }
  </style>
<body>

  <!-- implementation of the hosted D3- and TopoJson js-libraries -->
  <script src="http://d3js.org/d3.v3.min.js"></script>
  <script src="http://d3js.org/topojson.v0.min.js"></script>

  <!-- map creation -->
  <script>
  // canvas resolution
  var width = 1000,
      height = 600;

  // projection-settings for mercator
  var projection = d3.geo.mercator()
      // where to center the map in degrees
      .center([0, 50 ])
      // zoomlevel
      .scale(100)
      // map-rotation
      .rotate([0,0]);

  // defines "svg" as data type and "make canvas" command
  var svg = d3.select("body").append("svg")
      .attr("width", width)
      .attr("height", height);

  // defines "path" as return of geographic features
  var path = d3.geo.path()
      .projection(projection);

  // group the svg layers
  var g = svg.append("g");

  // load data and display the map on the canvas with country geometries
  d3.json("blok_zona.json", function(error, topology) {
      g.selectAll("path")
        .data(topojson.object(topology, topology.objects.blok_zona)
            .geometries)
      .enter()
        .append("path")
        .attr("d", path)
  });

  // zoom and pan functionality
  /*var zoom = d3.behavior.zoom()
      .on("zoom",function() {
          g.attr("transform","translate("+
              d3.event.translate.join(",")+")scale("+d3.event.scale+")");
          g.selectAll("path")
              .attr("d", path.projection(projection));
    });

  svg.call(zoom)*/

  </script>
</body>
</html>

and my JSON files: https://www.dropbox.com/s/xcp9iks661rzftx/blok_zona.json?dl=0


texts: Seeking algorithm to find point on surface with highest visibility?

texts:
I have a DEM which I want to find a point with largest amount of visibility ,I know that viewshed leads to visibility from a point but how can I find a point with highest visibility ?

One way is to loop every pixel in my DEM and then calculate visibility for each pixel. but this is not a very clever method.

Is there any tools or Algorithms which I can use to achieve this goal?

I can use c++ to run the algorithm.


texts: Downloading 1 arc-second elevation data for Canada other than Geobase?

texts:
Can someone suggest a source for Canadian elevation data better or equivalent to 1 arc-second?
I have a good data source for Ontario but for other provinces folks often recommend Geobase which claims that their elevation data ranges between 0.75 arc seconds to a maximum 3 arc seconds depending on latitude. However, I find that i cannot successfully generate shaded relief from the geobase data in ArcGIS. I've followed their instructions but the results are consistently un-usable. SRTM data from either cgiar.org or usgs.gov are good but I would like more detail.

Response from GeoBase:
Thank you for contacting GeoBase for a follow up on improvements to avoid 'stair stepping' effects found in certain CDED data sets.  Presently there is no change attributed to the status of CDED specification requirements.  CDED data originates from multiple partners and data sources which make it difficult to manage a homogeneous product throughout.
As you may now be aware the Z resolution in CDED files is causing the stair stepping effect.  According to CDED specification standards for this scale, all elevations readings are rounded to the nearest meter of elevation for the entire Canadian coverage with approximately an elevation point every 23 metres by 16-11metres for the 1:50K scale.   In flat areas, the rounding effect of the elevations causes homogeneous elevation plateaux instead of having a very gentle slope because the grid height values varies just a few centimetres.  It's these adjacent plateaux that cause the stair stepping effect.
In certain cases some smoothing effect may be achieved by using the Focal Statistics tool in ArcMap Spatial Analyst toolbox, neighbourhood tools before processing CDED files.  Depending on this software application you may be able to adjust viewing features with the tool menu for '' Neighbourhood=Circle, Radius=5, Units=Cell and Statistics Type=MEAN''.  But I don't think there is a perfect setting for all CDED files. There is a part of trial and error in the settings for each tile.  Because of this situation and the software you are using you may need to request assistance for your software provider to smooth out the look of you files.
Because no perfect solution was arised from our tests, GeoBase did not provide a universal user guide for smoothing out CDED elevation data, in the Frequently Asked Question segment.


texts:
Im a geomatics student,

During a geo referencing with QGIS, a .Points file is saved. I open it with Notepad++

It is written under this form :

mapX,mapY,pixelX,pixelY,enable
18.41330472103005533,-34.31455203862660852,18.41112659235668403,-34.30085589171973481,1

I would like to know what the different columns are?

The first coordinates mapX and mapY, are Daymark returned by the user.

Can anyone help me understand the meaning of these different columns?


texts: How can I draw isoline plots in ArcMap (ArcGIS Desktop) from various depths in water column?

texts:
We generally draw surface plots of variables such as nutrients and pigments from the surface of the water depth=0 but we also sample at other depths (10, 25, 30m) at different locations in the oceans.
If I have an excel file with 5 stations so 5 lat and long locations but each station has 5 depths and 5 corresponding nitrate concentrations -- How do I ask ArcMap to plot only that layer of nitrate between say 5-15 meters? We usually chose layer Properties->Symbology->Graduated symbols


texts: I'm unable to used drop_na() function for a data.frame object

texts:
library(dplyr)
library(tidyr)
library(mapview)

# species name Quercus leucotrichophora

gbif("Quercus", "leucotrichophora", download = T)

sp1 <- gbif("Quercus", "leucotrichophora", geo=T, download = T) #  download is essebtial for data.frame

class(sp1) # object is a data.frame

dim(sp1)

spg <- sp1%>%select( lon,lat)``` # extract only lat long columns

the  error it is throwing is as follows

Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘select’ for signature ‘"data.frame"


texts: Processing tool window doesn't appear in QGIS

texts:
I'm encountering a strange problem: I cannot execute any tools in QGIS (e.g., buffer, join--literally everything anything in the Vector or Raster menus. This seems to apply to all tools in all plug-ins, including the QGIS tools, GDAL, Grass, and Saga. I can see the tool options, but clicking one (e.g., Vector > Geoprocessing Tools > Buffer) deselects QGIS. It looks as if another window is opening but I am unable to see it or interact with it.

This error is happening 3.10.4, but also occurred in version 3.8 and 3.12. I have uninstalled and re-installed QGIS numerous times, and have deleted all leftover files in Program Files in between installs. The last build where QGIS worked as normal was version 3.8.

Any ideas?


texts: Making WMS GetFeatureInfo request gives Invalidpoint response?

texts:
When I make a GetFeatureInfo request to a WMS service, I get an invalid point response. I wonder what is missing from my request.

My Layer definition :

var Parcels = new TileLayer({
      name: 'Parcels',
      source: new TileWMS({
        projection: 'EPSG:4326',

        url: 'https://blabla.com',
        params: { 'LAYERS': ['myLayer'], 'TILED': true }

      })
    });

Map definition:

var mymap = new Map({
  layers: [Parcels],
  target: document.getElementById('map'),
  pixelRatio: 1,
  view: new View({

    center: fromLonLat([171.00, -45.00], 'EPSG:3857'),
    zoom: 9

  })
});

Map event:

map.on('singleclick', function (evt) {

          var layer = findByName('Parcels');
          console.log(layer);
          var url = layer
            .getSource()
            .getGetFeatureInfoUrl(
              evt.coordinate,
              map.getView().getResolution(),
              map.getView().getProjection(),
              {
                'INFO_FORMAT': 'text/xml',
                'propertyName': 'parcel'

              }
          );

When I click on a map it returns:

Invalid point (I parameter must be smaller than width parameter and J parameter must be smaller than height parameter.

It makes a request with coordinates higher than the Width and Height of the Tile map.

https://blabla.com?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=myLayer&INFO_FORMAT=text%2Fxml&propertyName=mylayerprop&I=706&J=282&WIDTH=516&HEIGHT=516&CRS=EPSG%3A4326&STYLES=&BBOX=-44.939169199785766%2C170.8956793139546%2C-44.93687999023354%2C170.8989132489774

I believe that this has to do with the openlayers getGetFeatureInfoUrl method options.

What am I missing here?


texts: Get WFS feature from click coordinates OpenLayers?

texts:
I have an OpenLayers map with a WMS background layer. When clicking the map, I need to draw the WFS feature located on the click.

Note that the vector layer is not loaded yet. I must load/request the mentioned feature on clicking.

How can I make a WFS getFeature() request with the click coordinates as a query parameter?

Is there is a better approach?


texts: Centroid latitude/longitudes of Landsat pixels

texts:
I was successful in getting centroid coordinates for pixels in Sentinel 2A image (10m by 10m).
But Landsat (30m by 30 m) posed a problem.

from osgeo import gdal
# Open tif file
ds = gdal.Open('C:\\Users\\ADMIN\\Desktop\\Data\\viveka.img')
# GDAL affine transform parameters, According to gdal documentation xoff/yoff
#are image left corner, a/e are pixel wight/height and b/d is rotation and is
# zero if image is north up.
xoff, a, b, yoff, d, e = ds.GetGeoTransform()

def pixel2coord(x, y):
"""Returns global coordinates from pixel x, y coords"""
 xp = a * x + b * y + xoff

 yp = d * x + e * y + yoff
return(xp, yp)

# get columns and rows of your image from gdalinfo
cols = ds.RasterXSize + 1
rows = ds.RasterYSize + 1
print(cols)
print(rows)
# Output and export
if __name__ == "__main__":
file=open("landsatfile2.csv","w") ## open file in write mode
for row in  range(0,rows):
    for col in  range(0,cols):
        x,y = pixel2coord(col,row)
        print "({} {})".format(x+15,y-15)
        file.write("{} {}\n".format(x+15,y-15)) ## write each coordinate separated

Why can't the code do it for Landsat?

The first image is from Sentinel and the second zoomed is from Landsat.
I want coordinates for centroid of each Landsat pixels.


texts:
This is my first post here so I apologize in advance if a convention is escaping me in the question format.

I am trying to display information from utfgrid in OpenLayers. I am using custom tiles in a XYZ format which I have setup as a baselayer. I have added a utfgrid layer using the openlayers demo but I cannot get it to work. The goal is to have information show up on mouse over events. The information from UTFGrid is passed to a div on the main html page.

Main html page:

main.js

var map;
    function init() { map = new OpenLayers.Map ("map", {

            restrictedExtent: new OpenLayers.Bounds(-157.3251,28.3044,-39.9032,76.5986).transform("EPSG:4326", "EPSG:900913"),
            units: 'm',
            numZoomLevels: 8,
            projection: new OpenLayers.Projection("EPSG:900913"),
            displayProjection: new OpenLayers.Projection("EPSG:4326")
        } );

        // This is the layer that uses the locally stored tiles
        var OLMCLayer = new OpenLayers.Layer.XYZ(
        "Local Tiles",
        "tiles/${z}/${x}/${y}.png",
        {
        type: 'png',
        isBaseLayer: true});
        map.addLayer(OLMCLayer);
        // This is the end of the layer

        //layer for interactive grid
        var utfgrid = new OpenLayers.Layer.UTFGrid({
        url: "tiles/${z}/${x}/${y}.json",
        utfgridResolution: 4, // default is 2
        displayInLayerSwitcher: false
        });
        map.addLayer(utfgrid);

        //callback
        var callback = function(infoLookup) {
            var msg = "";
            if (infoLookup) {
                var info;
                for (var idx in infoLookup) {
                // idx can be used to retrieve layer from map.layers[idx]
                    info = infoLookup[idx];
                    if (info && info.data) {

                        msg += " + info.data.csdname + ";

                    }
                }
            }
            document.getElementById("attrs").innerHTML = msg;
        };

        var utfControl = new OpenLayers.Control.UTFGrid({
        callback: callback,
        handlerMode: "move"
        });
        map.addControl(utfControl);

        map.addControl(new OpenLayers.Control.ScaleLine({'geodesic': true}));

        if( ! map.getCenter() ){
            var lonLat = new OpenLayers.LonLat(lon, lat).transform(new OpenLayers.Projection("EPSG:4326"), map.getProjectionObject());
            map.setCenter (lonLat, zoom);

        }

    }


texts: Weird streaks when working with a large raster

texts:
I am currently working with performing a large-scale landscape analysis in R and I am running into the following problem.
I have this large elevation raster that I made by using the mosaic function in the raster package. I use raster files from the ASTER database. Each tile represents a one-degree by one-degree portion of the earth with an estimated raster resolution of 30m^2. I include the code that I used to mosaic them below, but I do not believe that it is the problem since I have successfully used the same code before.
Mosaic <- do.call(mosaic,c(list of rasters, tolerance = 1, fun=mean)

Once that is run, I get the following elevation tif. I have included a picture of it and a link to download the file for the purpose of providing materials that aid reproducibility.

Link: https://drive.google.com/file/d/1yU_YDIqmCU-WqBUNftYiPPLpAzEMrQ8B/view?usp=sharing
Now, once I have the large raster file, I need to measure the topographic ruggedness of the raster and I use the terrain function in the raster package to do so. I use the Terrain Ruggedness Index as my method for measuring ruggedness
TRI <- terrain(Mosaic, opt="TRI")

This yields me the following raster.

As you can see in this example, there are faint streaks in the map which I cannot explain. Due to the full extent of my study, I need to aggregate this raster to a coarser resolution. Therefore I use the aggregate function in the raster package as well.
TRI_Agg <- aggregate(TRI,fact = 255, fun = mean)

This yields the following raster.
Here you can see that the streaks are much more pronounced and fairly consistent throughout the whole map. I tried doing this with other study areas and the problem persisted (i.e this would happen with raster tiles all around the world). I can't explain why those streaks are occurring. I have attempted different ways of going about this which include doing a neighborhood analysis (to fill in missing values if any) and processing the tiles before mosaicking them together. None of these have worked.


texts: Is the ArcMap Spatial Bookmark feature or something similar available in Inkscape?

texts:
When editing text in maps in Inkscape I wish I could have view/spatial based bookmarks as is available in ArcMap.

For example: Quickly zooming into the text at various parts of the map.


texts: Copying features to new location using QGIS?

texts:
I want to copy some features from one layer and paste them onto another layer with all their attributes, but in a different geographic location. Is this possible? These are very small features a very long way apart, so the cursor-based 'move feature' tool is not really an option.

I'm using QGIS 2.6.


texts: How can I determine the area of each block of a SpatialGridDataFrame in R?

texts:
I'd like to be able to determine the area of each block inside a grid of (long, lat) coordinates.  How can I do this?  Alternatively, how can I determine the area of any Spatial Polygon with a long+lat projection? (I can probably extrapolate to my grid scenario if someone can show me how to do this.)

Setup

library(rgdal)
library(sp)

# US bounding box
usbb <- c(left=-125, bottom=20, right=-60, top=50)

# Generate the grid vertices (Each vertex is 10 degrees away from the one above it and the one to the left of it.)
usgridDF <- expand.grid(Lon=seq(usbb[["left"]], usbb[["right"]], by=10), Lat=seq(usbb[["bottom"]], usbb[["top"]], by=10))
usgridDF[, "GridId"] <- seq_len(nrow(usgridDF))

# Convert from data.frame to SpatialPointsDataFrame
usgrid <- SpatialPointsDataFrame(coords=usgridDF[, c("Lon", "Lat")], data=usgridDF[, "GridId", drop=F])

# Convert from SpatialPointsDataFrame to SpatialPixelsDataFrame
gridded(usgrid) <- TRUE

# Convert from SpatialPixelsDataFrame to SpatialGridDataFrame
usgrid <- as(usgrid, "SpatialGridDataFrame")
proj4string(usgrid) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

To Do

# Determine the area in each block of the grid, in sq km
# ...


texts: mapbox with markers and checkboxes - disable all from showing on initial load

texts:
I have a page I built to display wind direction/velocity in various locations along a bike trail: http://microflush.org/cgi-bin/pathInfo.cgi

I have several markers in a legend, with check boxes to enable/disable them from showing up on the map.

My goal is to have NONE of the checkable markers show up on the initial load. However, when the page is loaded for the first time, ALL show up-- even though they're not checked. If I select 1, they all go away except for the one I selected, which will then disappear if I un-check that box.

Is there a way to make them all be invisible at first, and only show up if they're selected?

(You can view the source of the link above to see the code. I wasn't sure if I should post that all, since it's kinda long with all the points I've added.)


texts:
I have an empty field in every shapefile called "District" and I want to pull the filenames from each of those files into that "District" field.

Here is my code so far:

import arcpy
from arcpy import env
arcpy.env.overwriteOutput = "TRUE"
arcpy.env.workspace = "C:\Users\XXXXXX\Desktop\d2000"
worksp = arcpy.env.workspace
files = arcpy.ListFiles("*.shp")
for shapefile in files:
    arcpy.AddField_management(shapefile, "District", "TEXT")


texts: How to use maximum and minimum operator in raster calculator QGIS

texts:
I would like to normalize a raster, using the raster calculator in QGIS 3.x but I don't know how to use the maximum and minimum operator. I have tried the following sequence but it does not recognize it:

("Vuln_Erosion @ 1"-minimum ("Vuln_Erosion @ 1")) / (maximum ("Vuln_Erosion @ 1") - minimum ("Vuln_Erosion @ 1"))


texts: PostGIS: Querying a record table and explanatory table

texts:
Basically I have two tables one containing a collection of records and one table that contains extra info about features in the first. i.e

Records Table:
GID
SITE
SPECIES
GEOM

and the second table has
GID
SITE
SITE_desc

Essentially the first table has many records for the same site while the second table only has one record (in one row) for each site with a description.

What i'm wondering is if there is anyway i can query both tables (for a specific species for example) and have the field Site_desc be printed in a new column for each record where the sites are the same (SITE=SITE).

Is this possible?

Regards
Paul


texts: How to calculate nearest point to a linestring in PostGIS?

texts:
I'm following this question to calculate POI nearest a road(linestring). I'm able to calculate nearest points in the linestring but I'm not able to find the distance from the POI to nearest point(vertex) on the linestring.

Here's my code for finding the nearest POI.

CREATE TABLE road(id serial PRIMARY KEY, the_geog geography(LINESTRING,4326) );

CREATE TABLE poi(gid serial PRIMARY KEY, name varchar, city varchar, the_geog geography(POINT,4326) )

Values are:

INSERT INTO road (id, the_geog) VALUES (123, ST_GeographyFromText('SRID=4326;LINESTRING(85.280194 23.296728,85.281572 23.297479)') );

To calculate nearest point :

SELECT poi.name,poi.city,ST_AsTEXT(poi.the_geog), ST_Distance(road.the_geog, poi.the_geog)/1000.0 AS distance_km,ST_AsTEXT(road.the_geog::geometry) FROM road, poi WHERE road.id = 123 AND ST_DWithin(road.the_geog, poi.the_geog, 10000.0) ORDER BY ST_LineLocatePoint(road.the_geog::geometry, poi.the_geog::geometry), ST_Distance(road.the_geog, poi.the_geog);

If the linestring is represented by this: [85.280194 23.296728,85.281572 23.297479], I want the result like:

poi                   vertex                 distance_km
85.280194 23.296728   85.280001 23.299876    3
85.289673 23.291987   85.281572 23.297479    5


texts: the location of two points group with the same coordinates are not overlaped

texts:
I have two points data with the same coordinates, i add one of them by adding from Excel file and other one from CAD file, but in Arcmap the should be overlapped but it is not working. i selected the same projection system for them but still i can not solve this problem. please help me.
for eg. the green point is the point which i used GO TO XY and the blue one which is selected is the point selected from attributive table. i selected the same coordinate system for them but still the locations are so different.