texts: Correct way to use WGS84 (lat/long degrees) in QGIS

texts:
I created new project in QGIS (v3.2.1)
Added a OSM Standard layer with QuickMapServices plugin. Current CRS is EPSG:3857

Now I want to use map in degree units. As far as I understand I need to convert CRS to WGS 84. So I pressed the button:

Project Properties | CRS dialog box appeared. I filtered CRS by string '84'. There were two potential candidate geographic coordinate systems in that list: EPSG:4326 and IGNF:WGS84G

It looks like both of them do not fit because of the map distorted in both cases. The original view (EPSG:3857) was following:

After changing to CRS to EPSG:4326 it became look like:

How to solve my problem with degree coordinate systems?

The following is the image of layer context-menu:

It has not Save as element.


texts:
I am trying to convert a GRIB2 file to a GeoTIFF file using gdal.
My Commands:
gdal_translate -b 1 -a_srs EPSG:3857 some.grib2 some.grib2.tif
gdalwarp -t_srs '+proj=longlat +datum=WGS84 +ellps=WGS84 +units=m +no_defs' -te -130.103387 20.196426 -60.885910 52.807544 -tr 0.0174532925199433 0.0174532925199433 -r bilinear some.grib2.tif some.grib2.tif.reproj
I am receiving the following error:

ERROR 1: some.grib2.tif.reproj, band 1: Failed to compute statistics, no valid pixels found in sampling.

Source GRIB2 File (some.grib2): https://www.dropbox.com/s/zt1jzje64xvabbz/some.grib2?dl=0


texts: Any android application for field mapping with Google map in background

texts:
I need a android/iPhone app, through which I can draw point/line/polygon while survey only. The app should have facility to use Google satellite imagery or provision for uploading base maps.
Please suggest, if any.


texts: ArcGIS Data Reviewer - extract results using Python

texts:
I'm evaluating the ArcGIS Data Reviewer extension to see if it suits our needs.  I've created a Python script that will run some data reviewer checks and I can confirm that is successfully writes results of invalid values to the Reviewer Table.

From there, I can use the Data Reviewer 'Options' button in ArcGIS to manually export the table to a spreadsheet or to run a report.

I would like to be able to export this data from the reviewer table using Python to fully automate the process.  However, I cannot find any way of doing so.

Is there any way to extract the Reviewer Table data using Python?

(Note that I can see that there are a set of tables in the reviewer GDB that I could theoretically query, but they are a set of related tables that would take a lot of work to get the data out in the same format as the reviewer table - which I guess is really just a view over this set of tables.  Can I do it without trying to recreate the Reviewer Table from these other low-level tables?)


texts: SQL Server aggregate function in ogr2ogr command

texts:
I would like to define create a GeoJSON bounding box using SQL Server's EnvelopeAggregate function, embedded in an ogr2ogr command. I have the following command:

ogr2ogr -f GeoJSON envelope.json "MSSQL:server=SERVERNAME;database=DBNAME;trusted_connection=ye
s;" -sql "select TR, Geometry::EnvelopeAggregate(Shape) as Shape from gis.SampleTable where TR = '000280965
8' group by TR" "-s_srs" "+proj=tpeqd +lat_1=49 +lon_1=-110 +lat_2=49 +lon_2=-77 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no
_defs", "-t_srs", "epsg:4326"

The resulting GeoJSON file contains null geometry, despite the same SQL query working as expected in SQL Server Management Studio:

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },

"features": [
{ "type": "Feature", "properties": { "TR": "0002809658", "Shape":  "00000000010405000000009899992166D5408061325554EB21C100E00B936336E0408061..."     }, "geometry": null }
]
}


texts: How to save corrupt map service back to Service Definition (*.sd) file?

texts:
I have a corrupt ArcGIS Server Install and don't want reconfigure the existing services.

I understand that one can take a service definition file (.sd) and publish it as a map service on ArcGIS Server.

Is there a way to save a map service back to a .sd file? Basically Map Service TO SD File


texts: Identifying file format of web map layer?

texts:
I am interested in extracting data from the map layer available here.

The data looks like this (I inserted linebreaks for readability):

Malta_SAR normal #LineStyle000 highlight #LineStyle001 SAR Zones Lybia (NO
INFO IN IMO REPORT) 1 normal #LineStyle000 highlight #LineStyle001 SAR Zones
Lybia (NO INFO IN IMO REPORT) 1 0 #LineStyle00 1 11.526082,33.171136,0
11.531247,33.196851,0 11.537371,33.213983,0 11.566443,33.298784,0
11.585108,33.364755,0 12,33.916667,0 12.218734,34.175,0
12.40997642454452,34.37477266876127,0 23.583333,34.333333,0 24.1667,34,0
25.150691,31.669654,0 SAR Zones Lines normal #LineStyle000 highlight
...

Can anyone identify what file format this is, so that I can figure out how to read it properly?


texts: error setting up QGIS server with apache 2.4

texts:
I keep getting the 403 error and I have followed the directions from here.  Anyone know what's wrong?

Adding QGIS Server to an existing Windows web server

I set the httpd_gis.conf to require all granted

LoadModule fcgid_module modules/mod_fcgid.so

FcgidInitialEnv PATH "C:\OSGeo4W64\bin;C:\OSGeo4W64\apps\qgis\bin;C:\OSGeo4W64\apps\grass\grass-6.4.3\bin;C:\OSGeo4W64\apps\grass\grass-6.4.3\lib;C:\Windows\system32;C:\Windows;C:\Windows\System32\Wbem" DefaultInitEnv QGIS_PREFIX_PATH "C:\OSGeo4W64\apps\qgis" DefaultInitEnv QT_PLUGIN_PATH "C:\OSGeo4W64\apps\qgis\qtplugins;C:\OSGeo4W64\apps\qt4\plugins" DefaultInitEnv TEMP "C:\Users\mtenold\AppData\Local\Temp"

Alias /qgis/ C:\OSGeo4W64/apps/qgis/bin/

<Directory "C:\OSGeo4W64/apps/qgis/bin/">
    #SetHandler fcgid-script
    SetHandler cgi-script
    Options ExecCGI     Require all granted </Directory>

and I copied the httpd_gis.conf into the apache\httpd.d folder.  However when I try and run this request
 http://localhost/qgis/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities

I get the 403 error

If I type localhost into the browser, I see the screen that says "It works".  So I know the Apache server is working.


texts: Gauss-Krüger system to Unity cartesian Coordinate Conversion

texts:
I have to deal with coordinate conversation into my unity (a game engine which uses cartesian co-ordinate) project but I am zero in this field as it is involved math, geometry, and geography with lots of jargon to understand. Here is my case:
I am getting OBJ (object file) which is on Gauss-Kruger System. I am also getting this kind of numbers data (maybe Gauss Kruger Coordinate) with the 3d object file.
"crs" : "M34",
"parameter" : "5720"

"datum" : "21000101"

Per object data:
"x0" : "10393.00",
    "y0" : "336847.00",
    "z0" : "134.84",
    "dx" : "128.00",
    "dy" : "128.00",
    "dz" : "2.49",
    "x0original" : "11033.00",
    "y0original" : "336847.00",
    "z0original" : "147.44"

I want to know that how exactly I place the model as I have several objects so coordinate conversion should work between UCS(unity coordinate system) and Gauss Kruger System. Previously I am able to convert my unity co-ordinate to Lat Long via Unity GPS Conversion project.


texts: What version of SVG should I be saving my images with for QGIS? Any particular options I should use/avoid?

texts:
I have a illustration in Adobe Illustrator that I'm trying to export as an SVG file and add to my map layout using the Print Composer of QGIS 1.7.3.

When I add the image, I get strange artifacts that aren't in the Illustrator version. I tried a few different version of SVG in the Save As dialog, but they all end up generating the strange artifact.

What version of SVG should I be saving my images with? Any particular options I should use/avoid?

The orange block/white dots pattern below is the weird artifact:

Here are the Save As options for SVG in Illustrator CS 3:


texts: Trouble loading data into a newly installed QGIS

texts:
I have just installed QGIS 2.16 onto a laptop running Windows 10.

However, I'm having trouble loading any data into the program which crashes everytime I do this with the error message:

qgis-bin.exe has stopped working

Any ideas on why this is occurring?


texts: PyQgis show SVG Symbol on differnet value of specific attribute

texts:
I am working on a QGIS plugin. Main feature of this plugin will be changing the icons on a specific layer partly based on a specific attribute (it will be more exciting but you have to start somewhere). However, I succeed in changing all points at the layer to one specific SVG file.

After the ability of changing all points to one specific icon. I was looking around, how to change the code. But I am rather confused by the rule-based and categorized approaches. So which one should fit best to my problem and how to deal with different SVG icons. Because it looks to me that I only can use SVG icons by creating a specific layer and apply it to qgis.


texts: Display GeoTIFF single band image in Python

texts:
I have a single band GeoTIFF image which displays correctly in QGIS. However I tried viewing the same file using matplotlib and rasterio but I get a different output.
Is there something that I am doing incorrectly?
How do I display it like it is in QGIS?
I would also like to use a colour map for representing the different values of the single band with a color.
The EPSG of the GeoTIFF is 4326.
Here's my code:
import os
import rasterio
from rasterio.plot import show
from mpl_toolkits.basemap import Basemap
from matplotlib import pyplot as plt
from osgeo import gdal

if __name__ == "__main__":

    input_Dir = 'generated_46.tif'

    tiff_file = input_Dir
    raster = rasterio.open(tiff_file)
    print(raster.crs)

    bbox = raster.bounds

    show(raster.read(),cmap="gray")

Here's how it appears in qgis:

And this is my python output:

Here is the GeoTIFF description:


texts: Splitting layer with many polygons into independent shapefilles?

texts:
Is there any way to save many shapefile, which are not multipolygons, into shapefile, and then iterate over them?

I have shapefile with plot polygons and I want to have each one independent in order to build a grid for each one of them and to use rotate in order to make the grid look good inside the polygon.

This is what I have right now

My end goal is to use the Model builder to iterate those polygons and to create for each one of them a grid.


texts: Error when using SAGA's "topographic correction" algorithm in QGIS

texts:
I'm trying to do topographic correction of Landsat images using SAGA's "topographic correction" algorithm from QGIS' processing toolbox. I was able to do it using Minnaert Correction with slope and cosine correction, but when I try C-correction I get the following error:

ERROR 1: Can't load requested DLL: C:\OSGEO4~1\apps\saga\dll\gdal_MrSID.dll
127: Impossibile trovare la procedura specificata.

and no image is created.
Anyone has any idea about how to solve this?

In alternative I wanted to use GRASS' "i.topo.corr", but it seems like it's not available in the Processing toolbox?

I am using QGIS 2.10.


texts: Using azimuths and bearings in GeographicLib?

texts:
I am currently using the JavaScript port of the geodesic routines from GeographicLib:

https://geographiclib.sourceforge.io/html/js/

https://github.com/yurijmikhalevich/node-geographiclib

I am computing the inverse problem where given lat1,lon1,lat2,lon2 the various azimuth's and the distance is computed:

var geod = GeographicLib.Geodesic.WGS84, r;
r = geod.Inverse(lat_aid, lon_aid, lat_loc, lon_loc);

The description of the function is given here:
https://geographiclib.sourceforge.io/html/classGeographicLib_1_1Geodesic.html

My use case is one where one set of my lat_loc,lon_loc is a given location, and the other set of lat_aid,lon_aid is a navigation aid used to get to the location. I am looking to calculate the distance and bearing from the aid to the location.

Unfortunately I am not a GIS expert or very versed in geography in general and I am have a few questions:

I am having trouble figuring out which azimuth I require (azi2 vs azi1)
Are the calculated azimuths based on grid or true north? Is there any real difference in this model since it estimates the Earth's shape.
The values I get seem to range from 0 to 180 and 0 to -180, am I correct in assuming that a positive azimuth value is clockwise from north, and a negative value is counterclockwise from north? ie: azi -70 would be a bearing of N 70 W and an azimuth of 70 would be N 70 E


texts: How to create Bookmarks in the Print Composer?

texts:
I try to find a way to create a bookmarks in the QGIS Print Composer, so when i move the map with the pan option i will quickly return to the wanted frame of the Print Composer.

I search in the view menu in the Print Composer but didn't found any solution.


texts: Is there a way to tell geoserver to send cp1250 http header?

texts:
I'm running geoserver with some web pages in data_dir/www. These pages
are encoded in cp1250, but when I load them using geoserver
(localhost:8080/geoserver/www) they are encoded or they behave as they
were encoded in utf-8. As a result i get tons of question marks
instead of real letters.

I am not able to convert pages into utf-8, so is there a way to tell
geoserver not to send utf-8 http header and send cp1250 header
instead?


texts: Basemap widget to access map layers from ArcGIS Server?

texts: Can a basemap widget retrieve a list of map layers from an ArcGIS server instead of ArcGIS Online or Portal? If not, how can I get or build a similar widget?


texts: Counting frequency of specific value in selected columns of an attribute table of the feature class in ArcGIS 10.4

texts:
I want to calculate the number of members in the houses (feature class which contain information of the members in the house interms of age,sex,occupation) based on the information entered in age column (if the column contain -9999 then there is no person in the house) as shown in following table

For this task i have used Model Builder as shown below but due to the NULL values generated in the process of calculation, I could not get the exact output.
Are there any ways to do this using Python?


texts: Can I move points by editing x/y in attribute table?

texts:
Im new to QGIS and I just cannot figure out how to do this one hopefully simple thing:

I have a bunch of points in shapefile and i need to move some of them around by hand / coordinates. I've understood, that if i move points by hand its possible to update their new xy to attributes table by field calculators $x and $y operators (Just select the column and insert that operator i suppose?).

However, how about reverse? I thought that the features position on screen could be changed by editing the xy columns on attribute table...but the point remains on the same old location. Am i doing something wrong?

thanks!


texts: Auto incrementing every second row in ArcGIS Desktop?

texts:
How do I auto increment every second row (field type is short and I work with 150 fields) in ArcGIS Desktop 10.3.1?

I saw this code in
https://support.esri.com/en/technical-article/000011137

rec=0
def autoIncrement():
 global rec
 pStart = 1
 pInterval = 1
 if (rec == 0):
  rec = pStart
 else:
  rec += pInterval
 return rec

 autoIncrement()

and the unwanted result is the blue arrow, whereas i need the result will be the as the yellow arrow (i changed it manually):


texts: Elevate points from imported DXF file to heights in LAS databse

texts: Is it possible to elevate points that are on the same plane in dxf to heights that are in LAS? I have imported dxf and las file in Arcmap but I need guidance or at least direction where to start looking for solution.


texts:
I'm using GeoServer since couple of months and it were working fine. Today, after some gap, when I typed localhost:8080/geoserver; it is giving an error

Object not found!

The requested URL was not found on this server. If you entered the URL manually please check your spelling and try again.

If you think this is a server error, please contact the webmaster.

Error 404

localhost
20/12/2013 5:51:07 PM
Apache/2.2.17 (Win32) mod_ssl/2.2.17 OpenSSL/0.9.8o PHP/5.3.4 mod_perl/2.0.4 Perl/v5.10.1


texts: Altering order of attributes in table using ArcPy

texts:
I have gone through a sequence of table joins. The resultant table should be with the attributes rearranged in more logical order. I cannot find this functionality within ArcGIS or any mention of it in Python

For example, say a TableResultingFromJoins has the following attributes:

v1,v2,v3,w1,x1

... and I want it like this:

v1,w1,x1,v2,v3

What do I do?


texts: Simulating flood without flooding whole area and covering surface structures (with DEM)?

texts:
I have a high resolution DEM and want to do a flood simulation starting at a river. My hope is that the water runs into historical (now dry) river channels (in which I'm interested in) at first and step by step. In my imagination the whole old channel system appears at the end.
But the problem is, the higher the water raises the area which is covered by water gets bigger and bigger. Like in this pic (on the right):

But I would love to get a result like on the left side.
Now I hope i can describe this understandable: My idea is to do a flooding step by step and every step will be "saved". At the end you have several flood phases which you can intersect.

Red, green and yellow are the different flood phases and the brighter blue is the intersection.
I dont know how to manage this (even after internet research) or if it's realizable at all. Maybe you have an idea or even a better problem solution!? I use ArcGIS with extensions.


texts: Is there a reason for using PostGIS functions such as ST_Area and ST_Length when there are already functions/operators for doing the same thing in PostgreSQL (area and @-@)


texts:
I have been working with GIS for many years but I am fairly new to the spatial analyst side of things. I am working with many satellite tracks for many different individual animals and would like to determine a migratory corridor for these animals. We have decided to use line density for the analysis of these data. Our question is: how to determine the search radius, which I have learned is very important?

I found out that one way to do this would be to use the tools on movebank.org like ArcMET and Home Range Tools. However to use these tools, a time/date stamp is required for every point on the satellite track and we do not have this information for all of the points. This is also the reason why we decided to use line density instead of kernel density. So we decided to use line density in ArcGIS.

I found that one way that many people determine a search radius for point density or kernel density is by using the Incremental Spatial Autocorrelation. However to perform this type of spatial statistics, one needs count data/incidence data. Since these are satellite tracks, the data inherently is only one animal per track (line).

So based on these issues, does anyone have any ideas on how to determine the search radius? Any other spatial statistics that we could use?


texts: MapInfo incorrectly counting lines intersecting a buffer

texts:
I made a buffer for all objects on start and end nodes of polylines then disaggregated them. I then used update table to derive the number of lines intersecting:

However this appears to be incorrect as a result of 8 is returned for the following buffer, when there are only 6 polylines:

Where am I going wrong here?


texts: QGIS expression using Aggregate and start_point to select a vector within a buffer

texts:
I am trying to put a default value of a field so that the value always updates once any change is make item. Either the expression builder or function editor are fine.  I'm more familiar with the expression builder so I went that path. There should only ever be 1 pole within 5m of the start point so I can use concatenate but if there is a way to find the closest pole that would be best.

If I use the code below I get no info no matter how large I make the distance.

 aggregate( layer:='Poles', aggregate:=' concatenate' , expression:="Pole Id", filter:=buffer(start_point( geometry(@parent)) ,0.00005),',')

With this code I get every item no matter how small I make the distance.

aggregate( layer:='Poles', aggregate:=' concatenate' , expression:="Pole Id", filter:=geom_to_wkt(buffer(start_point( geometry(@parent)) ,0.00005)),',')

I tried to find an answer but this was the closest I could find.  I don't have enough rep to add a comment so I had to post a new question. @csk  Close but not quite answer


texts: Reducing the processing time when handling millions of 3D points to 3D lines in PostGIS

texts:
I am trying to calculate the minimum distance between 16 millions of 3D points and 3D lines.

The query is working for more or less 1 million data, but after that, it took so long and never ends. All the tables have spatial index.

The query is:

SELECT
    public.points.orden,
    st_length(st_shortestline(public.points.point_geometry, public.line01.geom)) as line01_length,
    st_length(st_shortestline(public.points.point_geometry, public.line02.geom)) as line02_length,
    st_length(st_shortestline(public.points.point_geometry, public.line03.geom)) as line03_length,
    st_length(st_shortestline(public.points.point_geometry, public.line04.geom)) as line04_length
ST_3DLENGTH(ST_3DSHORTESTLINE(public.points.point_geometry, public.line01.geom)) as line01_3dlength,
ST_3DLENGTH(ST_3DSHORTESTLINE(public.points.point_geometry, public.line02.geom)) as line02_3dlength,
ST_3DLENGTH(ST_3DSHORTESTLINE(public.points.point_geometry, public.line03.geom)) as line03_3dlength,
ST_3DLENGTH(ST_3DSHORTESTLINE(public.points.point_geometry, public.line04.geom)) as line04_3dlength
into public.lengthdist
FROM public.points, public.line01, public.line02, public.line03, public.line04
LIMIT (SELECT COUNT(*) FROM public.points);

After, I select the min distance:

SELECT public.lengthdist.orden,
    min(line01_length),
    min(line02_length),
    min(line03_length),
    min(line04_length),
    min(line01_3dlength),
    min(line02_3dlength),
    min(line03_3dlength),
    min(line04_3dlength)
FROM public.lengthdist
GROUP BY public.lengthdist.orden

And also y set the endpoint for each min line with ST_EndPoint, and I have the same problem.

I tried with CREATE TABLE - INSERT INTO statement as well, same results. Also with an UPDATE. Also, I tried with just ST_SHORTESTLINE without ST_LENGTH.

How can I do to reduce the never-ending processing time?

I am working in Windows 10 x64, 16GB RAM.

DB Version: 11
OS Type: windows
DB Type: desktop
Total Memory (RAM):16 GB
CPUs num: 4
Connections num: 20 Data Storage: ssd
max_connections = 20
shared_buffers = 512MB
effective_cache_size = 4GB
maintenance_work_mem = 1GB
checkpoint_completion_target = 0.5
wal_buffers = 16MB
default_statistics_target = 100
random_page_cost = 1.1
work_mem = 22573kB
min_wal_size = 100MB
max_wal_size = 2GB
max_worker_processes = 4
max_parallel_workers_per_gather = 2
max_parallel_workers = 4
max_parallel_maintenance_workers = 2

P.D.: I've VACUUM ANALYZE all tables too.


texts: Create network from points in QGIS using grouping variables

texts:
There seem to be some answers for this problem for ArcMap, but not for QGIS.

I have a set of points representing the power grid and I want to produce a network from them. Here's my data:

ID     lat_from    lat_to   lon_from    lon_to    voltage
1        51.1       51.2      0.2        0.97      400
1        51.2       51.9      0.96       0.84      400
1        51.4       51.8      0.5        0.37      275
1        51.6       51.7      0.7        0.57      400
1        51.7       51.6      0.6        0.67      275
1        51.8       51.4      0.9        0.97      400

Many of the tools don't seem to accept directional connections (from point x to point y) or grouping variables (e.g. voltage)

I'm looking for something like this:

Are there any suggestions for how to do this?


texts:
I'm part of a group making global visualization software. I need to find earth imagery (visual range and hopefully IR as well) between 500 m/pixel and 75 m/pixel.

The blue marble next generation data doesn't get high enough resolution.

There seem to be some pay services for this like Planet Observer, TruEarth , and Digital Globe. I'd like to avoid paying if I could.

I see that there is Sentinel 2 data which I think might meet my needs but I'm having trouble finding how to download it. Even on the ESO website I can't find a  good way to do it en masse (or at least more than one image at a time.

In the end I need the data tiled so that my project can zoom. I could probably do that myself if the data is an a reasonable form to begin with.


texts: Creating TOC widget for RTL languages in ArcGIS API for JavaScript?

texts:
I'm trying to make the TOC widget to run from right to left in ArcGIS API for JavaScript.

Is it possible to achieve, and how?


texts: Opening and saving multiple shapefiles as .tab files?

texts:
How do I open several shapefiles, save as Tab files then append as one table using MapBasic code?

I have created a list that contains the tables and path to tables.  But cannot loop to save tables individually.

I have reviewed many of your posted answers (really good answers) but still cannot find solution.


texts: Calculating and outputting max and average distance from point to points and output to table?

texts:
I have two point layers. One contains nest locations. The other contains movement points for the individuals being tracked. Both tables have a individual id parameter (one nest per id, many movement points per id).

I want to create a script (ModelBuilder preferred) that will generate a table (or update the nest table) to show the maximum distance and average distance between an individuals nest and the movement points from the other table.

So far I've tried using the PointDistance tool to create a table of distances from the nest to the points. Then I've had to manually get the max and average. As I have to rerun this often, I'd like to have the script automatically get the max and average and update the table or create a new table with one row per nest ID showing the max and averages.

The first image is a basic table layout and what I'd like in the output. There is really no relevant data in the tables other than the ID for the individual. I'm also including an image showing a sample nest point (in green) with the movement points for that nest ID.

I am using ArcGIS 10.4.1. License type: Advanced.


texts:
I tried by following the instructions described here but they don't work. It looks like the repository mentioned in those instructions is actually not working:

  http://fedora.vitu.ch/QGIS/7/x86_64/repodata/repomd.xml: [Errno 14]
  HTTP Error 404 - Not Found

I was unable to find any other working method. This other repo I was able to find only goes up to CentOS 6: http://elgis.argeo.org/


texts: Merging overlapping GeoJSON LineStrings (Tippecanoe) and add visualization

texts:
I would like to merge overlapping GeoJSON LineStrings and visually style merged LineString. I have a database of GeoJSON LineStrings which represent routes people drove. I would like to visualize routes which people drive most frequent with different colour than the ones that are not that frequently driven.

Example of Tippecanoe rendered Mbtiles of GeoJSON LineStrings (numbers stand for overlapping routes):

And I would like to achieve something similar to what I depicted on image bellow:

I came across few articles how to visualize GeoJSON Points, but I can not find any helpful tip how to achieve the visualization for LineStrings.

Technologies that I'm using: MongoDB, Tippecanoe (https://github.com/mapbox/tippecanoe), tileserver-gl (https://github.com/maptiler/tileserver-gl)

Ideally I would like to achieve similar results that Strava has:

I don't know if technologies I'm using are good fit to achieve that.


texts: Displaying each polygon segment while digitizing polygon with Polygon construction tool?

texts:
I am almost positive that I've digitized the polygons and that I have seen each segment shown as I dropped in vertices during the sketch.  I figure that I inadvertently configured something and now I can't get it back.

Once I finish the sketch (by double-click or F2), the polygon redraws and displays as expected.  The symbol has an outline and there is only one feature template.

Here's the polygon, mid-sketch:

I'm expecting lines for the first 2 segments, but not getting them.

Here's the polygon after I've finished the sketch:

Here is my configuration while digitizing the polygon:

And this is the construction tool I'm using:


texts: QGIS run python script upon print composer export

texts:
I'm looking for a way to execute a python script upon exporting a print composer. Specifically, I want it to run once for each page exported. (logging some data and possibly updating a database).

So far I've created a custom expression that writes a log to a file and returns an empty string. I created an off-page label that runs the expression.

The problem is, this runs every time the page refreshes so it's running a bunch just paging through the atlas.

Trying to get it to run only on export, I put it in the 'export filename' field but that too runs every time the page refreshes.

Is there some way I can either force it to only run on export, or in the python code detect if the page is being exported before running?

edit:
here's the test code I'm running so far:

def logTest(logitem):
    global logFileName
    f = open(logFileName, 'a')
    f.write("{}\n".format(logitem))
    f.close()

@qgsfunction(args='auto', group='Territory')
def postData(postItem, feature, parent):
    logTest(str(postItem))
    return u''

and my expression to test it is simply:

postData("id")


texts: PostGIS ST_LineMerge: why do 3 lines not merge, but any 2 do?

texts:
In PostGIS, ST_LineMerge is used to "stitch" several lines into a single linestring. If it cannot, it returns the original multilinestring.

I have a scenario where I have three lines in a set. They happen to be over a river. When I attempt to use ST_LineMerge(geom[]) with three lines, it does not merge. However, using any two of the three lines does work!

Here is some code:

select ST_asText(ST_LineMerge(
  ST_Union(
    Array[
      -- THIS IS THE EAST-MOST LINE.
      -- The remaining lines start where this line ends.
      ST_GeomFromEWKT(
        'LINESTRING(-110.597626872104 45.0142399967956,-110.6000500721 45.0151951967941)'
      ),

      -- THIS IS THE MIDDLE SHORT LINE.
      -- NOTE: It falls just short of the third line.
      ST_GeomFromEWKT(
        'LINESTRING(-110.6000500721 45.0151951967941,-110.602446129524 45.0154430716088)'
      ),

      -- THIS IS THE LONGEST LINE.
      -- It appears to "contain" the second line.
      ST_GeomFromEWKT(
          'LINESTRING(-110.6000500721 45.0151951967941,-110.602478272096 45.0154463967937,-110.604730272093 45.0151811967941)'
      )
    ]
  )
))

Note that the second line is likely what's causing this problem. Each of these lines are (assumed to be) subsections of a river linestring. I believe that the middle line would be considered by PostGIS to be a subsection of the third line, even though they are not sharing more than one point. Indeed, the second line falls just short of the middle of the third line. When I replace the second coordinate of the second line to share the second coordinate of the third, it merges correctly.

Success in this case means these three lines are merged into one, and the new line's length is equal to the long western linestring plus the short east-most line. In other words, the new linestring is the simplest, shortest line possible.


texts: Clear/refresh table of contents in Geoprocessing tools

texts:
I'm building a Python toolbox and when I run one of my tools, I noticed that the table of contents doesn't update itself in the drop down list of my Python tool. As you can see in the image below, some layers are present in the list and are not in the table of contents although I've updated my data frame (by the way all of these layers are created by a previously executed Python tool).

These are the functions I used to update my data frame, without any result:

arcpy.RefreshActiveView()
 arcpy.RefreshTOC()

Also, the increment after the GPLx layers adds up every time I re-execute the tool. This means that somewhere in the system, an identical layer is stored. That's why Python increments the layer number suffix (I guess).

I thought that ArcGIS probably stores the layer history somewhere so I looked for it but didn't find it.

EDIT

Run the following tool and put any layer you have in the input parameter. Once it's done, remove the layer from your data frame and go back to re-execute the tool. You will notice in the drop down list that the removed layer is still there:

# -*- coding: utf-8 -*-
import arcpy

class Toolbox(object):

def __init__(self):
    """Define the toolbox (the name of the toolbox is the name of the
    .pyt file)."""
    self.label = "LayerTest"
    self.alias = ""`

    # List of tool classes associated with this toolbox
    self.tools = [LayerTestTool]

class LayerTestTool(object):

def __init__(self):
    """Define the tool (tool name is the name of the class)."""
    self.label = "LayerTestTool"
    self.description = ""
    self.canRunInBackground = False

def getParameterInfo(self):
    """Define parameter definitions"""
    params = []

    param_in = arcpy.Parameter(
            displayName="In Layer",
            name="Layer",
            datatype="GPFeatureLayer",
            parameterType="Required",
            direction="Input")
    params.append(param_in)

    param_out = arcpy.Parameter(
            displayName="out",
            name="out",
            datatype="GPFeatureLayer",
            parameterType="Derived",
            direction="Output")
    params.append(param_out)

    return params

def isLicensed(self):
    """Set whether tool is licensed to execute."""
    return True

def updateParameters(self, parameters):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""
    return

def updateMessages(self, parameters):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""
    return

def execute(self, parameters, messages):
    """The source code of the tool."""
    Layer = parameters[0].valueAsText
    out = arcpy.MakeFeatureLayer_management(Layer, "myoutput")
    arcpy.SetParameterAsText(1, out)
    return


texts:
I need to convert a pair of (latitude, longitude) coordinates to a POI if available. I am able to reverse-geocode the coordinates to an address. After that is there a way to map the address to a POI or anything else ?

Could someone point me to free API if available ?


texts: QGIS Convert vector classification map into raster

texts:
I tried to convert vector classification into raster file for accuracy analysis following these suggestions

My vector data is like this:

However, none works
The first one, interpolation, resulting in .asc file which is invisible with zero value.

The second one I tried to use rasterize function, it loaded, but nothing happen, the output is not change. I tried to follow this video tutorial by copying style, but the raster layer become invisible.


texts: Deleting input vector or raster being used in QGIS 3 modeler

texts:
As it is explained in the title I am writing a script to process this task.

The script is already written but I need it to be triggered when the inputs are not needed anymore. In other words : I need to set the parents algorithms box in my script. Like in the illustration below :

Illustration comes from this post.

It seems for the "processing.run" method we can set that with the is_child_algorithm=TRUE parameter?

But I use "pure" Python in my script so it doesn't matches.

I can't find the method for my case in QGIS documentation.

Does anyone has the corresponding code for this case ?


texts: Running ERDAS Imagine and ArcGIS Desktop on same machine?

texts:
ArcGIS Desktop 10 and ERDAS 8.7 were already installed in my computer. But I want to open both programs at the same time. If ArcGIS can open properly, ERDAS cann't open and vice verse, showed that following error (Fig.1).

IF I switch off ERDAS using License Manager Tools "LMTOOLS" (Fig.2), then ArcGIS can open properly, If I switch off ArcGIS using the same way above, ERDAS can open properly.

So what I want to ask are what is the reason for this problem and how to fix it?

I've changed "any" into 27002 and both "this_host" and "son" into "user_PC" which is my computer full name. This time Erdas can open correctly, but arcgis cannot. The following image is a problem with opeining ArcGIS Desktop.

Can you tell me what this means?

Thank you for your answer, U2ros. Here I'm adding some more photos of license files. According to your advice, when I opened both license files for arcgis and erdas, they look like this. First image is ArcGIS license file which is very similar with you described above except "son". I think I can change 27001 into 27005 etc.

Next image is license file of Erdas. This looks totally different from you described above. Here I'm fully confused which one is port number that I need to change to something like 27002. Can you give me further advice to for this.

.


texts: QGIS 3.4 handle bad layers and PostGIS views

texts:
My project contains several views from PostGIS (2.3). I had to update them and now the Handle Bad Layers window shows for just three of them.
I checked the data source and the query after table=... is the same in PostGIS. I have just added a primary key in the views with row_number() over(). I removed it and load the project again without success.
By the way, I noticed that in this last case (removing pk), the data source in the project is not updated so I edited the .qgs to do it. There it is the same stored in the project and in PostGIS. It seems to have no effect.

The layers have many styles and the project many layouts so you can imagine the pain if I can't fix the data source issue.


texts: QGIS Hidden attributes not hidden in Relations

texts:
I have a PostGIS view, which a referencing layer for another layer for many to one relation in QGIS. I have hidden the following fields id, table_name, oid of the postgis view so that the user would not see those fields in the attribute table. It does work in form view. The fields are hidden as in the image.

But they are visible in the attribute view. I tried Hide Column option in the attribute table but it did not work. Is there any way I can hide them both in form view and attribute view?


texts: updating the extent automatically in the print composer

texts:
Just started using QGIS (after trying what I need to in GRASS).

Using the Print Composer, I've put down a map which is based on the map in my map canvas, which is what I'll want to print. Then put down a few logos as well, which need to be printed on the map.

Now, the people who will eventually be using this map will be panning and zooming around merrily on the map canvas, then expect it to print what is on the map canvas, as well as the logos. Using a saved Print Composer file here will be fine.

Trouble is, can the Print Composer automatically update it's extents to match the map canvas? If not, I don't think it's a HUGE problem - I'll just tell them to use the "Set to map canvas extent" button, but it would be nice if it could do this automatically.


texts: Trying to get PRECISE vector coordinates using geom.asPolgon

texts:
I've got a bit of a problem. A script I've written is trying to get the coordinates of every vector in a layer, which is MAINLY working:

output_file = open('c:/myfile.txt', 'w')
layer = iface.activeLayer()
for f in layer.getFeatures():
  geom = f.geometry()
  if geom.type() == 1:
    if geom.isMultipart():
      line = '%s, %s, %s\n' % (f['USRN'], f['Street_Nam'], geom.asMultiPolyline())
    else:
      line = '%s, %s, %s\n' % (f['USRN'], f['Street_Nam'], geom.asPolyline())
  elif geom.type() == 2:
    if geom.isMultipart():
      line = '%s, %s, %s\n' % (f['USRN'], f['Street_Nam'], geom.asMultiPolygon())
    else:
      line = '%s, %s, %s\n' % (f['USRN'], f['Street_Nam'], geom.asPolygon())
  unicode_line = line.encode('utf-8')
  output_file.write(unicode_line)

However, there seems to be an issue with the Y coordinates (in my case, I'm using WGS 84 (EPSG:4326). The longitude is coming out okay, but the latitude is not precise at all. For example:

-0.0920572,51.6096

Obviously, the longitude is fine, but the latitude isn't. This is an excerpt from a shape which is an oval, and the longitudes are all different. However, because of the lack of precision, the latitudes that are output for this feature are all the same (51.6096) which, when I put this into my software, renders it as a single line.

Sorry, forgot to mention I'm using QGIS 2.8.1.


texts: postgis table - shapefile - mapinfo file

texts:
I have a .tab file with a ridiculous amount of columns (96). One of my colleagues at my new job wanted to update them manually one by one. I told him I could do that in an easier way using PostGIS. So I told him to give me the data and let me do my thing.

I converted the .tab file to a shapefile, uploaded the shapefile to PostGIS, uploaded a table to PostGIS with the data which was needed for the update. I updated all the 96 columns at once and achieved the wanted ending in a table in PostGIS. He wanted it back in a .tab file.

I exported the file as a shapefile, which went perfect. Then I export the file as a .tab file using QGIS. So I figured I had the end result and have send it to my colleague. He checks it and reports back to me that the last 10 columns are empty, which they should not be. Also the coordinates differ and the polygons are not in the same location as they were in the original file in MapInfo, which he uses.

I checked it in QGIS and there is no problem at all. The columns are filled, the polygons overlap with the original data. I am totally perplex.

I know that my coordinate system RD new sometimes differs in code. I remember adding a small line of text in the formula to align it once. So I guess the movement of the coordinates is something I can fix myself. But why do columns appear empty while they are not in QGIS?

And why can some .tab files I export in QGIS not be opened in MapInfo? It gives an error the moment I open the .tab file.

I am researching myself, but I need some guidance to find the right direction to research further to. I'd love some help and that is why I ask the question here.


texts:
Is there any backup strategy for GIS data?

Say .mxd, shape files, data base connection files, feature layers, geo-database etc.

According to ideal methodology:

1.On external hard drive : best to take backup of all folders and sub folders.

2.Working on ArcGIS JSAPI so taking daily backup on CVS

Similar fashion what would be the recommended backup strategy for GIS Data?

We are using ESRI and open source platform (Desktop and server based + spatial database)

Any help would be great!!


texts: How to write Python if/then statement using information from two fields in field calculator

texts:
I am using ArcGIS Desktop 10.6.1. I create a new field of type short integer in the Attribute Table. My new field is named CIPOT and I want to populate it with values from the field CI, but only if values in field Z = 100, otherwise I want CIPOT to equal zero.

My Python script:

def NumValue( CI ):

  if (Z == 100):
    return CI
  else:  return 0

NumValue(!CI! )

Here is a screen shot:

I am getting an error message that global name 'Z' in not defined.


texts: Unite several polygons that share boundaries in QGIS

texts:
My study zone is in the middle of two geologic maps with two different shp files.

Is it possible to remove the boundaries of the polygons between these two maps that are created at the union of these in order to have only one polygon (per geologic unit)?.

In other words: to obtain an "unpartitioned" map (in QGIS).


texts: Copy & paste between layers crashes program

texts:
I attempted to copy a feature in one layer and paste it into another layer, as I have often done. However, I got the message, "Crash dumped; minidump written to c:\Users..." The program then closed itself.

I uninstalled the program, restarted the computer, and reinstalled the latest version, 2.14.21.

Upon starting the program, I get the yellow exclamation mark message, "Plugin processing: The plugin will be disabled because it crashed QGIS during last startup...re-enable the plugin when the problem has been solved."

When I click on "Plugins" in the main menu, I get a list, and the explanation that ones showing in red have a problem. However, none of mine are red.

When I again try to copy and paste from one layer to another, I still have the problem—even using different layers.

Any idea what's going on, and what I might do about it? I am running Windows 10.


texts:
Is there an ArcMap equivalent to the Sieve Tool found in QGIS?

  Sieve Tool:  This utility removes raster polygons smaller than a
  provided threshold size (in pixels) and replaces them with the pixel
  value of the largest neighbor polygon.


texts: Unexpected result using ModelBuilder in batch?

texts:
I am running a model which I created with ModelBuilder at 10.3. The model is supposed to make a selection by attribute on a feature layer, and use that layer to clip a raster. The model is working fine when I run it, however if I try to use it in batch mode it gives unexpected results: I have many files with the same extent and content as the first line of the loop (of the batch).

The parameters of my models are:

the input feature layer against which a selection is specified
the input raster (.BIL)to be clipped with the selected feature
the SQL expression for the selection (of the type "FIELD" = 'value')
the output folder and filename (with .BIL extension)

Basically, in the batch window the feature and raster parameters are fixed, I am only changing the SQL expression to point a different 'value', so that the output raster is clipped on different areas each time, and the filename name of the output.

Below is my model:


texts: Error when trying to plot legend for geom_point with ggplot2

texts:
Here's some mock code to go along with my problem:

#read in libraries
library(ggmap)
library(sp)
library(raster)
library(ggplot2)

#read in country shapefile, project, fortify for mapping with ggplot
tz <- getData('GADM', country = "TZA", level = 1)
tz.prj <- spTransform(tz, CRS("+proj=longlat +datum=WGS84 +no_defs"))
tz <- fortify(tz.prj)

#create dummy data frame
lon <- c(36.68299, 36.66667, 31.80501, 33.87462, 35.75161,
         32.23135, 39.20833, 36.68173, 39.26093)
lat <- c(-3.386925, -3.216667, -1.329641, -2.018532, -6.162959,
         -2.885038, -6.792354, -8.139894, -6.822502)
tz.c <- data.frame(lon, lat)

#plot
ggplot() + geom_polygon(data = tz.prj,
                        aes(long, lat, group = group),
                        fill = "#f1f4c7",
                        color = "#afb38a") +
  geom_point(data = tz.c,
             aes(lon, lat),
             pch  = 21,
             color = "black",
             bg = "red") +
  coord_equal()

Output looks like this:

I want my legend to include the red dots from geom_point to define what they are (let's just say hospitals for this example). While my country outline is a solid color right now, it will eventually be symbolized in some other way with the hospital points on top of it, so I'd like to include all legend items in my legend.

I did some searching and found that in order to use show.legend = TRUE with geom_point, you have to plot everything within aes(). When I try to do this, it gives me the error: Error: A continuous variable can not be mapped to shape.

I converted my lat and lon to factors to try to quell this issue to no avail. So, how do I work with this data so that I can plot it to look the same as the image I've provided while also being able to plot/define the legend?


texts: Open .NET implementation of the OGC Geometry Object Model from Simple Feature Access 1 specification?

texts:
I am looking for a .NET implementation of the geometry object model described in OGS's Simple Feature Access, Part 1 specification.

Preferably lightweight. I don't need WKT support, just the capability of instantiating spatially referenced geometries and the ability to determine basic spatial relationships between them.

Any suggestions?


texts: QGIS QTiles / QMetaTiles poor image quality

texts:
While making slippy map tile sets from a QGIS (2.10.1) project running on Vista 32, I've noticed some reduced image quality from QTiles 1.5.0 and QMetaTiles 1.0, versus the output from a rendered print composer (and subsequent fairly intricate workflow to end up with tiles).

In all cases, I'm exporting to .png, 256x256; the checkbox for 'Make lines appear less jagged...' was on for all these test cases, and the image quality setting (with .png selected, this is only adjustable in QMetaTiles, not QTiles) makes no discernable difference.

A few images are attached.  The aliasing on the text is one thing, but the stairstepping and staggering on the contour lines is the main issue.

Can anyone explain why these differences occur, and/or point out a way to get the 'better' version using the QTiles or QMetaTiles plugins?

Composer workflow:

QTiles:

QMetaTiles:


texts: VB Script field calculator - error with conditional statements

texts:
I have limited VB Script experience, but have been working on the script below and keep getting errors, the most common being "expected identifier".

Can someone tell me what I'm missing here?

Dim class as Integer
if [AGE_DATABA] >= 200 AND [Riparian_S] <> 1 then
class = 5
elseif [AGE_DATABA] > 200 AND [Riparian_S] = 1 then
class = 2
else
class = 0
End if

expression = class


texts: Counting unique features of points inside polygon in QGIS

texts:
I have a series of points with various fields in the attribute table, including a car_id. I also have several polygons in which some of these points are located.

I want to count how many unique car_id's occur in each polygon. For example, if there are multiple points in a polygon where car_id = 25, the car_id = 25 should be counted once. The car_id count for a polygon using the sample data below should be 4.

car_id = 25
car_id = 25
car_id = 25
car_id = 22
car_id = 22
car_id = 23
car_id = 24

What is the best approach to achieve this? Or is there an appropriate QGIS tool available in QGIS 3 to do this?

I know that in QGIS 2.x the tool 'Count Unique Points In Polygon' would have worked. However this tool is not available anymore in QGIS 3.10.


texts: Create PostGIS-Topology (edges, nodes, faces) from a list of edges

texts:
I have a large in memory table of edges with the following structure:

facenr1, facenr2, x1, y1, x2, y2

x1,y1 to x2,y2 is a line.
facenr1 and facenr2 is the facenumber of the left and right face.

Now I want to load this into a PostGIS-Topology scheme, so that I create a correct topology and eventually get the faces in the faces table.

The documentation states that I can create this using INSERTs:

http://trac.osgeo.org/postgis/wiki/UsersWikiPostgisTopology#LoadingTopologydata

How do I do this?


texts: How can I connect a Microsoft Access Database and a FGDB, such that information put into the former is automatically transferred to the latter?

texts:
I am working on a project for the National Park Service. The goal of the project is to connect a Microsoft Access Database to a File Geodatabase, where any information put into the .accdb tables will automatically be put into the FGDB tables. The reason for this is that it will enable other employees in our park without GIS experience to easily put in data from cave surveys into the simplified UI we made for the accdb. We want that data to then be transferred to the FGDB as seamlessly as possible; the end goal is to have something that gets it done with one process and no other hassles involved.

So far I have connected the access database to ArcCatalog via an OLE Database link, which works without issue. However, I am having trouble figuring out how to set up a way in which data from the accdb will be automatically sent to the FGDB, without having duplicates that we would have to manually delete. The only way I could think of was to set up a model like in image I included below; it only includes the model for just one of the table datasets needed.  It almost works, but it keeps deleting the wrong set of duplicates. It deletes the new set (the most current data), not the old set. Essentially, I am just trying to have the model overwrite the older cave survey data and replace it with the newer data. If I could find a way to have the model delete the older duplicate set, it would do exactly what I want it to.

Can anyone give me input on my method or recommend another?

I have limited python experience.


texts: Using arcpy.GetParametersAsText() into Text file?

texts:
I've just started using Python for ArcMap. I'm trying out a script that works before I put arcpy.GetParameterAsText(0) into it. Does anybody know what I'm doing wrong? (I put the parameters on the command line).

Coding is for Icelandic letters

#-*-coding: cp1252 -*-

import arcpy
from arcpy import env
import os

k = arcpy.GetParameterAsText(0)
arcpy.env.workspace = k
inFeatureClass = arcpy.GetParameterAsText(1)

txtFile = open(os.path.join(env, "GGSkrár15.txt"), "w")
txtFile.write("Upplýsingar um skrár fyrir gagnagrunn Póst- og Fjarskiptastofnunar" + "\n")
txtFile.write("-------------------------------------" + "\n")

desc = arcpy.Describe(inFeatureClass)
desc2 =  "Skráarheiti: {} Shape: {} Gerð: {}\n".format(desc.name, desc.shapeType,
                                            desc.datasetType)
txtFile.write(desc2)

desc3 = arcpy.Describe(inFeatureClass)
desc4 = "Extent: {} Hnitakerfi: {}\n".format(desc3.extent, desc3.spatialReference.name)
txtFile.write(desc4)

field_list = arcpy.ListFields(inFeatureClass)

for field in field_list:
    line = "Dálkaheiti: {}, Gerð: {}, Fjöldi stafa leifður: {}\n".format(
         field.name, field.type, field.length)
    txtFile.write(line)

txtFile.close()
print "Script completed"


texts: Fast Clustering Algorithm for Geographic Data?

texts:
I have been browsing through the net in my attempt to find a fast clustering algorithm to cluster geographic points, but lot many of the results found were for clustering algorithm in general and not for geographic data.

Some of these algorithm expected me to specify the number of cluster in advance which is something I don't want to do?

Any suggestion to explore some algorithm?


texts:
I have fetch data from database and draw it on qgis 2.8 as vector layer it works fine with limited data but for large data it will not work and goes into not responding stage, I am unable to understand the problem. My code is follows

 layer =  QgsVectorLayer("Point?field=gid:integer&field=lattitude:double&
 field=longitude:double&field=depth:integer
 &crs=EPSG:4326","items","memory")

 QgsMapLayerRegistry.instance().addMapLayer(layer)

 for i in range(len(Receivers)):
      # add a feature
      feature = QgsFeature()
      X = Receivers[i][1]
      Y = Receivers[i][2]
      temp=[]
      #X= X cordinate  longitude
      #Y= Y cordinate  latitude

     feature.setGeometry(QgsGeometry.fromPoint(QgsPoint(Y,X)))

     length=len(Receivers[i])
     for j in range(0,length):
           temp.append(Receivers[i][j])
     feature.setAttributes(temp)
     layer.startEditing()
     layer.addFeature(feature, True)

     layer.commitChanges()
     layer.updateExtents()

Receiver list contain the data fetch from database (postgres)


texts: Getting just land and water layers in QGIS with QuickOSM?

texts:
How can I just get land and water layers that I can style with different fills and outlines?

Really - Just "land is green, water is blue" layers (that I can then add other layers on top of)

I'm wanting/tempted to use QuickOSM and get a vector layer of one or many polygons to use as "land" such that I can color the land and also color the water (or background, or "not land").

I'm new to GIS so please let me know if there are better ways to do this. Natural Earth Isn't detailed enough - I'm looking at roughly 1:10000 scale.

I've been able to get some boundaries and lines and admin zones but none are just the land.  They all are bits and pieces that would need to get stacked and flattened and there must be an easier way (right?).  I'm exporting as an image into illustrator/inkscape so I'm looking to create areas that I can fill and style and outline etc.


texts: Parsing and filtering JSON to CSV using ogr2ogr flags

texts:
I want to export filtered data from JSON to a CSV by selecting a couple of fields and applying a where clause. I get empty output when using either the -sql or -where parameters with the -select (field_list) parameter.

Here are two different formulations of the command:

ex 1)

ogr2ogr -f CSV "output.csv" "raw.json" -select "id","name","type" -where 'type="city"

ex 2)

ogr2ogr -f CSV "output.csv" "raw.json" -sql 'select id,name,type from raw where type="city"'

The following query successfully exports desired fields but unfiltered.

ogr2ogr -f CSV "output.csv" "raw.json" -select "d","name","type"


texts: Listing Feature Classes of Multiple geodatabase in Multiple Folder?

texts:
I'm new in python scripting and try to list polyline feature classes of multiple Geodatabases. The Geodatabases are in many nested folders.For example in folder "A" I have 20 folders and in each folder i have 3 geadatabases. the question is that how can i list the feature classes of the geodatabases using Arcpy. I know simple Python codes to list feature classes in a folder such as :

import arcpy
from arcpy import env
env.workspace = ".././Data/6"
fclist = arcpy.ListFeatureClasses()
print fclist

I surf in this site and found a  Q&A about How to list feature classes in multiple Geodatabases in a folder? . It was very useful but only for multi Geodatabases not multi folders.


texts: Buffer multiple roads to find points within and compare?

texts:
What I am trying to do is buffer a main st and a side street and get all of the points that intersect that buffer, however they need to be mutually exclusive. In the example below, my results would give me a duplicate at the intersection of Epping Forest Rd and Point No Point Dr.

Logic: I need to buffer Epping Forst Rd and get all points that intersect that buffer and buffer Point No Point Dr and get all points that intersect that buffer then compare both results and drop any points that are duplicate.

So far I have this SQL statement:

SELECT  aa_points.int_id
    FROM    aa_points,
            road_ln
WHERE road_ln.st_nam = 'EPPING FOREST'
AND ST_Intersects(aa_points.the_geom_webmercator,
ST_Buffer(road_ln.the_geom_webmercator,10))

That queries one road but I am at a loss on how to do the second road in a nested statement and how to compare/drop records.


texts: Downloading static high resolution images using Google Earth

texts:
I wanted to download some high resolution (<5m) images. I have their co-ordinates and viewed them on https://earth.google.com/web/ and it's exactly what I want. I have approximately 300 such images to download. I was wondering, how can I do this with Python?
I checked out the Google Earth API but I am not sure which dataset to choose from the catalog to get the same images from the web version. Any idea as to how I can do this?


texts: how to extract overlap from non-dissolved buffer zones in qgis 1.7.4

texts:
I have a shapefile with the inventory of all limetree in a region. I made a buffer around each tree.  So, the place where buffers overlap each other is the best place to put beehives.

Is there a way to extract the overlap between non-dissolved buffers in the same layer?  Or is there another way?

Thanks for the help.
Pieter


texts: How to fit Google Earth satellite image to UK OS 50 DEM

texts:
I have extracted a DEM tile from the UK Ordnance Survey's free OS Terrain 50 dataset
The tile contains spot heights for a grid comprising 51 nodes (West to East) by 51 nodes (North to South) with inter-node spacing 50m giving areal coverage of 2500m x 2500m.

I can load this DEM data into a home-made, browser-based, 3D viewer application (built with WebGL, THREE.js, javascript).

Now I wish to copy (for private use only) a satellite image from Google Earth (non-pro version) and apply it as a texture to the 3D viewer app.  I can do this in a rough and ready way by screen-snapping an image from Google Earth, then pasting the image into a paint program and then cropping and re-sizing.  I have done this manually in a rough way but getting accurate location and boundaries for the image will be laborious.  Also I expect that there may be horizontal positioning errors due to differences in the projections used by OS Grid and Google Earth.

I would like to know if there is a better way of extracting and processing the satellite image (using free software) to fit it to the DEM grid?


texts: Openlayers custom (local) OSM tiles layer with limited extent

texts:
I know, how to add an overlay layer by specifying my own URL. But what if the overlay layer extent is not as large as base layer - e.g. overlay layer is only avaliable for one country and not for the whole world.

  new OpenLayers.Layer.OSM("My custom local overlay layer",
  ["http://a.tile.myserver.org/maps/${z}/${x}/${y}.png",
  "http://b.tile.myserver.org/maps/${z}/${x}/${y}.png",
  "http://c.tile.myserver.org/maps/${z}/${x}/${y}.png"]);

How can I pass parameter, that OL would not seek for the tiles, that does not exist - that are outside of the avaliable layer extent (outside that country)?


texts: How to extract records with a specific but shared attribute

texts:
My goal is to isolate the ID#s observed in Zone C and only Zone C (like IDs 2236 and 3007 in the attached table). There are records for which the same ID was observed in other zones besides Zone C, or was observed multiple times in Zone C.

How can I extract only those individuals who were observed in C and ONLY in C? Selecting by location is not what I need because some of the ID#s in C are also in other zones.  I have tried Select by Attribute and the Dissolve tool but I am not quite getting the right combination. Can someone point me in the right direction?  I am using ArcGIS Desktop 10.4.


texts: Is it possible to export raster data from PostGIS or MapServer as a base64 string?

texts:
Is it possible to output a PNG in a base64 string so that can be directly loaded into a browser? I've read through the PostGIS documentation which refer to the GDAL output formats. The GDAL output formats here state the following:

  GDAL includes support for reading, and creating .png files. Greyscale,
  pseudo-colored, Paletted, RGB and RGBA PNG files are supported as well
  as precisions of eight and sixteen bits per sample.

Regarding MapServer, the only reference I can find is an unanswered question asked on the mailing list six years ago.

I am wondering if I can generate a png text string in the following form so that it can be loaded into a browser using the Data URI scheme:

<img alt="" src="data:image/x-png;base64,iVBORw0KGgoAAAANSUhEUgAAABAAAAAQCAYAAAAf8/9hAAAACXBIWXMAAAsTAAALEwEAmpwYAAAABGdBTUEAALGOfPtRkwAAACBjSFJNAAB6JQAAgIMAAPn/AACA6QAAdTAAAOpgAAA6mAAAF2+SX8VGAAABLElEQVR42qSTQStFURSFP7f3XygyoAwoYSYMPCIpk2egMFSmUvwCRpSRDIwYGbwyVuYykB9y914m951z7nHe6J26dc9u77XXWmdvJLF7/audqx9JYuvyW92LL0li8K2df2r17CPEVk7ftXTclyQqAMmRCwC5I3fS42a4W7y74VYDNAAuJA8AaXIsSACsDgAdAJeFrnnyoMBygKZJJ3b1It0AmsTMDPdEgrujJqHEwCxqznMaD2KgyCDRnEuo8qJhHvx/hcQDbzGoix5Yi4G1TcwZWNEDKwJU+WDkhg2ToDaD+M65YcVB8jg3Y5IY5VQAyyf9gLJw+CqAuYNnAczsPQpgevtBU937kDexcdssj8Ti0ZskMd97CRs3u//U2sjJzbtwH1+/Cf8jS/gbAMmWc42HzdIjAAAAAElFTkSuQmCC" />

I would like to test this approach, as I am interested in the approach that MySociety used when creating Mapumental. In the description of how they made this tool, they refer to color-cycling of a base64 PNG string. The approach used by Mapumental is very complex - I am just interested in exploring how to even start this approach.


texts: MapServer WMS with inline Element export is not exported correctly in print composer

texts:
I have a WMS service from MapServer 7.0.4 with a trademark as POINT inline Feature in the bottom corner defined with image coordinates. Therefore, it is always in the same place.

As soon as I export the map with print composer, due to the tiling, depending on the resolution, the trademark is displayed several times.

The problem is present from QGIS 2.1.4 - 3.0.2.

Has someone an idea how to solve this?

This is the MapServer code for the Trademark

This is the prind composer window.

This is the exported Raster image (600dpi)


texts: Preventing ArcPy from adding new features to Current Workspace?

texts:
The code below takes a huge raster, and a set of point features. For each of them creates a circular buffer (100), then squares it, and then clips that part of the buffer. In doing so, we end up with a png image for each of the point locations.

The script works fine, the main problem, is that it adds each clipped feature to the workspace. These are around 35,000 points, so at about 1100, it just kills arcmap, and I have to restart the process from where it stopped (that is the reason there is an if, to check for 1107).

Is there a way to just do the clipping, without everything being loaded in the environment?

import arcpy
import os
image_name = "MyMosaicDataset2.tif" #has to be loaded in the env
main_path = 'C:\Users\leon\Documents\Data\VCR_Data\Leon'
workspace_name = 'Mars_Tartarus_Colles_Georefrencedg.mxd'
arcpy.env.workspace = os.path.join(main_path, workspace_name)
base_for_buffer = 'VRCs'
buff_name = 'VRC_buff'
output_buffer = os.path.join(main_path, buff_name+'.shp')
output_sq_buffer = os.path.join(main_path, buff_name+'_sq'+'.shp')
radius_buffer = '100 meters'
arcpy.Buffer_analysis(base_for_buffer, output_buffer , radius_buffer, 'FULL', 'ROUND', 'NONE','','PLANAR') #This is to create the buffers
arcpy.FeatureEnvelopeToPolygon_management(output_buffer, output_sq_buffer,"SINGLEPART") #This transforms the circular buffers into rectangular buffers
rows = arcpy.SearchCursor(buff_name+'_sq')#here are the buffer polygons, and stored ina  list
#This is for clipping
cnt=1
for row in rows: #go over the list
    if cnt >= 1107:
        feat = row.Shape
        name = 'vrcp_'+str(cnt)
        output_raster = os.path.join(main_path, name)
        arcpy.Clip_management(image_name, "#", output_raster+'.png',feat, "0", "ClippingGeometry")
    cnt+=1
    enter code here


texts:
Has anyone played around with using ArcObjects in Python at version 10.1?

I am having all sorts of problems, such as TypeErrors when wrapping the ESRI OLB modules and AttributeErrors when attempting to use a method or property of a wrapped module.

Anyone had any better luck?


texts:
I have a code that will read a .tif and split the image into a specific number of rows and columns (code 1). I modified it to be able to be able to easily set the wanted image size (by pixels) (code 2). The wanted size wont be perfect because I am using the celi function to round to the nearest int. That is not an issue. Basically I am modifying QGIS way of creating a grid and clipping to each cell.
With the split images, I am running various classification methods to detect specific things. Unfortunately I cannot share my data with you. The problem that I am facing is that the left side of my images are not being split and therefor no images are able to run with the classification. You can see how in image 1 I have my classification areas and in image 2 the classification areas do not reach all the way over. There is a specific section where they stop. And yes, there are areas here that would be identified.
I've tried this method on a few images now and the same thing occurs. It's always on the left side. I have set my rows/columns to 3 and 3 and the missing section is no longer missing. It does split. So I do not understand why when I increase my columns and rows to a image size of say 600x600 (pixels), I end up missing data.
Is there something in my code that I am missing or does anyone have any suggestions?
# Setting the directory
os.chdir(r"D:\ortho")

# Loading in the image
rimg = gdal.Open("image.tif")

# Upper left corner of the minX and maxY
gt = rimg.GetGeoTransform()

xmin = gt[0]
ymax = gt[3]

res = gt[1]
xlen = res * rimg.RasterXSize # units is important UTM or WGS
ylen = res * rimg.RasterYSize

# how many tiles you want to have in each row and column
xdiv = 70
ydiv = 105

# Determining the size of each new image
xsize = xlen/xdiv
ysize = ylen/ydiv

print(xsize)
print(ysize)

xsteps = [xmin + xsize * i for i in range(xdiv+1)] # plut because we start in the left top corner where X is at its lowest
ysteps = [ymax - ysize * i for i in range(ydiv+1)] # minus because we start in the left top corner where Y is at its highest

for i in range(xdiv):
    for j in range(ydiv):
        xmin = xsteps[i]
        xmax = xsteps[i+1]
        ymax = ysteps[j]
        ymin = ysteps[j+1]

        # Splices the image up into the set divs and saves them.
        gdal.Warp("D:/data/image" + str(i)+str(j) + ".tif",rimg,
                  outputBounds = (xmin,ymin,xmax,ymax),dstNodata = -9999)

wanted_x = 2000
wanted_y = 2000
y_pixel_len = wanted_y * res
x_pixel_len = wanted_x * res

# how many tiles you want to have in each row and column
xdiv = math.ceil(xlen/x_pixel_len)
ydiv = math.ceil(ylen/y_pixel_len)

# Determining the size of each new image
xsize = math.ceil(xlen/xdiv)
ysize = math.ceil(ylen/ydiv)

image 1

image 2

Edits
There are some warnings that occur when I run the code. I think these are just images that contain no data. So I end up getting a black image. I do have code that will delete those images. Not sure if this helps???
Warning 1: for band 1, destination nodata value has been clamped to 0, the original value being out of range.
More than 1000 errors or warnings have been reported. No more will be reported from now.


texts: Rasterio applies a random error to all pixels when saving a raster

texts:
When I read a raster and then save it with rasterio the pixels are sligthly modified. Code example:

import rasterio

inputFile = './band3.jp2'
with rasterio.open(inputFile) as b3:
    v3 = b3.read(1)
    print('Before: ',v3)
    profile3 = b3.profile

##Saving raster
output='./test.jp2'
with rasterio.open(output, 'w', **profile3) as dst:
    dst.write(v3, 1)

##Reading the raster I just saved
with rasterio.open(output) as b3:
    v3 = b3.read(1)
    print('After: ',v3)

This outputs:

Before:
         [[  0   0   0 ... 665 632 661]
         [  0   0   0 ... 728 687 704]
         [  0   0   0 ... 703 702 700]
         ...
         [  0   0   0 ... 613 605 603]
         [  0   0   0 ... 610 600 605]
         [  0   0   0 ... 608 606 607]]

Warning 6: driver JP2OpenJPEG does not support creation option TILED #I don't think this is important

After:
        [[ 11   0   0 ... 663 629 655]
        [  0   0   0 ... 717 680 687]
        [  0   0   0 ... 689 693 687]
        ...
        [  1   0   0 ... 604 597 587]
        [  3   0   0 ... 604 595 593]
        [  0   0   0 ... 594 594 589]]

I'm not sure if it is important because I didn't modify it but the profile looks like:

{'driver': 'JP2OpenJPEG', 'dtype': 'uint16', 'nodata': None, 'width': 10980, 'height': 10980, 'count': 1, 'crs': CRS.from_epsg(32615), 'transform': Affine(10.0, 0.0, 600000.0,
   0.0, -10.0, 1600020.0), 'blockxsize': 1024, 'blockysize': 1024, 'tiled': True}


texts: Replicating the "Calculate Geometry" table tool in python

texts:
I'm trying to automate a very tedious process in our company that involves switching the data frame's coordinate system and performing calculate geometry on a couple of fields, these are the two situations I'm having trouble with :

First situation would be !SHAPE.CENTROID.X.! using data source's coordinate system with the unit conversion specified in the picture

Second is switching the data frame coordinate system, choosing to use it and perform the following calculation :

So, two main questions :

1) How to switch the data frame's coordinate system and calculate based on it

2) Get !SHAPE.CENTROID.X.! with the specified unit on picture 1

I'm using arcpy.CalculateField_management for this since there already is a selection

I really searched for this one but couldn't find what I'm looking for


texts: Is the result of "generalization" algorithm the same every time using it in ArcGIS?

texts:
I have a database with roads, and i update it oftenly. For another department we give them the same database but generalized so they can fit the forest to the roads.

I wanted to know if you use generalize on the same database several times, do you get the same results regarding lines that haven't change since?

I mean, do the unchanged roads have still the same geometry and vertex number so they still fits with the forest from the previous work?


texts: ArcObjects snipped to Select all Features currently visible on the ArcMap screen (ActiveView)

texts:
I guess I have to get the extent of the ActiveView and then use it to perform a spatial search.

Is that the proper way to do it?

Please share the code snipped if you have any.


texts: Rasterio showing a blue (blank) image wrong

texts:
I am trying to open a raster in my Jupyter Notebook using rasterio. I began with two Landsat 8 images, which I processed in QGIS to calculate the NDVI. I want to plot the original images used for the calculation (band 4 and 5) along with the NDVI.
My NDVI plot however is a blue image, in QGIS and Windows however it does show correctly.


texts: Find single nearest point of point in mysql

texts:
I have an app that formerly used PostgreSQL with PostGIS but now the environment changes and I have to use MySql.

In my app i had a database query that searches the nearest point to a point. Note that my table has about 40 Million rows.

SELECT
    strength,
    created,
    updated,
    point,
    ST_Distance_Sphere(ST_GeomFromText('Point(:lat :long)'), point) AS distance
FROM signalstrengths
ORDER BY distance ASC
LIMIT 1;`;

I read that mysql can't use postgres' <-> operator and the above query is what I came up with through reading docs and SO.

But this query doesn't work and I'm concerned that the speed of the query will suffer from changing to MySql (from what I've read).

Can you help me fixing that query and tell me more about the speed differences?


texts: I am using “Hindi and Gujrati” font then I am getting “Black Squre Dots”?

texts:
At the time of generating PDF file when I am using “Hindi and Gujrati” font then I am getting “Black Squre Dots”. Below I am using code. Here I am trying to convert in UTF-8 but it’s not converting.

file_path_PDF_new="C:\example\hello.pdf"
doc = SimpleDocTemplate(file_path_PDF_new,pagesize=letter, rightMargin=72,leftMargin=72, topMargin=72,bottomMargin=18, encoding='utf-8')
Notice = "नमस्ते दुनिया"
new_Notice=Notice.decode('utf-8')
styles=getSampleStyleSheet()
styles.add(ParagraphStyle(name='Justify', alignment=TA_JUSTIFY))
ptext = '%s' % new_Notice

Story.append(Paragraph(ptext, styles["Normal"]))
Story.append(Spacer(1, 12))

doc.build(Story)


texts: "TypeError: list indices must be integers or slices, not tuple" - the code works in one case but not in another?

texts:
The following code executes smoothly:

for i in range (1, 10):
    for j in range (1, 10):
        if a1[i][j] == 100:
            if np.count_nonzero(a1[i-1:i+2,j-1:j+2]) > 8:
                a1[i][j] = 1

While the following gives the mentioned error:

for i, item in enumerate(myarray):
    for j, element in enumerate(item):
        if na2[i][j] == 1:
            if numpy.count_nonzero(na2[i-1:i+2,j-1:j+2]) > 2:
                na3[i][j] = 0

I am not able to understand the reason.


texts: Exporting map layout to spot color EPS file in ArcGIS Desktop?

texts: File export option allows to export to cmyk eps file. I need to export a map layout to spot color EPS file in arcgis 10.2.2


texts: Questions about protected areas' effectiveness assessment

texts:
I am researching about the protected areas' effectiveness assessment,
My question is for example one protected area is built in year 2001, with MOD13Q1 data available since year 2000, how to better translate the  effectiveness of 'protected area'? I know one of the method can be done that is the comparative analysis of 'before protected' and 'after protected', like: 2000-2001----before protected, 2001-2019: after protected, but for 'before protected' the length of time is just one year and this result is not objective enough. I want to know is there any good solutions to resolve this problem of the extremely short former period: 'before protected', is there any way or mainstream method to substitute this unreliable method?


texts: Setting character encoding for all layers at once using PyQGIS?

texts:
I have to set the character encoding for all layers from system to UTF-8.

Is there a way to change this in all loaded layers at once?

I'm using QGIS 2.8.1 on Windows 7.


texts: How to delete label style saved in ArcView?

texts:
The answer to my question may be very simple, but I can't seem to find it myself with online searches

I have symbology and text/label styles (etc) that I've saved for easy access when I'm working on a Map.  This is mostly for reasons of consistency.    If I wanted to delete one, how would I go about accomplishing this?

EXAMPLE (see image below):  Highlighted in yellow, are duplicate label styles I saved at some point.  I want to do some house cleaning.  How do I delete one of them (for dummies please)?


texts:
I am using Proj4js to transform between WGS84 and a local Swedish grid RT90 (EPSG:2400). It works in FF, but in IE the Proj objects often is not readyToUse when I need them to.

I use this code to initialize the transformation objects (part of a bigger object):

wp.js

tera.wp = {
reproject: {

    referenceSystem2400: null,
    referenceSystem4326: null,

    initializeReferenceSystems: function() {

        Proj4js.defs["EPSG:2400"] = "+proj=tmerc +lat_0=0 +lon_0=15.80827777777778 +k=1 +x_0=1500000 +y_0=0 +ellps=bessel +units=m +no_defs";

        var self = this;
        if (self.referenceSystem2400 == null) {
            self.referenceSystem2400 = new Proj4js.Proj("EPSG:2400");
        }
        if (self.referenceSystem4326 == null) {
            self.referenceSystem4326 = new Proj4js.Proj("WGS84");
        }

        var loaded2400 = self.referenceSystem2400.readyToUse;
        var loaded4326 = self.referenceSystem4326.readyToUse;

        alert("init: " + loaded2400 + ", " + loaded4326);

    }
},

init: function() {
    var self = this;

    self.reproject.initializeReferenceSystems();

    // In this call, the transformation is performed
    this.repopulateFromSession(function() {

        self.disableLoader();
    });
},

the web page

<script src="http://maps.google.se/maps/api/js?v=3&sensor=false&language=sv&region=SE" type="text/javascript"></script>
<script src="/scripts/proj4js/lib/proj4js.js" type="text/javascript"></script>
<script src="/scripts/proj4js/lib/defs/EPSG2400.js" type="text/javascript"></script>
<script src="/scripts/wp/wp.js" type="text/javascript"></script>

<script type="text/javascript">

    $(document).ready(function() {
            tera.wp.init();
    });

The explicit declaration of the 2400 parameters is to really ensure that the parameters are ok, and does not have to be loaded from anywhere. (I've been at this a while now)

Although, when I load the page and run the initializeReferenceSystems function, the alert shows init: false, true indicating that the 2400 object has not loaded properly.

What can I do to ensure that both my transformation objects are ready to use like this?


texts:
I am very new to ArcGIS Desktop.

I have a jpeg showing an old map of NYC zones. I was able to georeference this image with a shapefile of NYC blocks. I wish to now create a new feature class where I can trace these zones from the old jpeg onto the new feature class layer using the blocks shapefile as a reference. I will also wish to trace along the blocks in that blocks shapefile layer.

I've been following this thread's advice (second reply) so far. It is as follows:

  Add your empty feature class/shapefile to arcmap, and click your Editor button and go to Start Editing. (Assuming you are using V10?)

  If you plan on following existing geographic features, look at using the Trace tool.
  Hint: When tracing, whatever feature you want your line to follow, it has to be Selected.

  There are several ways to select a feature:
  ...by Location
  ...by Attribute
  ...by Graphic

  Try to explore the Feature Construction toolbar and the various segment creation tools.

However, after I create this blank feature class and add it as a layer, I cannot edit it! I go into "Start Editing" mode, choose the blank feature class, but I cannot trace or do anything after that.

Is there any hope for me?


texts:
Using ArcGIS Basic 10.1, there is a symbol of a paper scroll next to the tool, instead of the hammer symbol.  That is the only tool I click on that crashes my program.

I am not sure what that paper scroll symbol means next to it.  Is there an update needed somewhere?  How else can I do the Point to Polygon conversion that I do need (without using scripting).


texts: Add new column in raster attribute table

texts: I would like to add a new column in a raster attribute table and populate it with values based on another column using if statement, the new column should have letters such as A, B, C and D. How do I do that in R?


texts: Creating global Geoprocessor object for ArcObjects application?

texts:
I have an application that does many operations including some Geoprocessor  tasks.
Should I create a single Geoprocessor() object for the lifetime of the application or should I create a new one each time I want to use a Geoprocessor tool?


texts: QGIS script accessing old version of layer

texts:
I have added a couple of fields to a postgis table which is loaded into a layer in QGIS. When I access the attribute table in QGIS I can see the two new fields but when I run a script (using scriptrunner) to access the layer it does not show the two new fields.

I check the fields in the script by listing them:

print [field.name() for field in layer.pendingFields() ]

but the new fields are not there and I get key errors when I try and access the fields by name or index.

How can a script access something different to the QGIS gui?

QGIS 2.18.7 on a Mac.


texts: Insert a mixed geometry type directory of shapefiles in a workflow with FME

texts:
I can very easily deal whole TAB or DWG files directories (for example) using the "Single Merged Feature Type" option in the respective Reader. It takes all geometries encountered and allows me to use the GeometryFilter to deal with each geometry separately after reading files.

  "Allowed geometries : All" is displayed in the reader

Surprisingly, I can't do this with shapefiles as easily. Every time I give FME a whole directory with different types of geometry shapefiles, FME takes the first feature he encounters and ignores the others or generates empty features.

  "Allowed geometries : shape_pointz" is displayed in the reader (the first layer encountered was points). I guess it does "Schema from first feature '*.shp [SHAPE]'"

The only way to deal with my respectives layers is to check "Individual Feature Types" but it generates as many readers as different layers, which is not at all what i'd like to do (the same workflow as described here by Mark Ireland) : FME Workbench: Merging Multiple Shapefiles into One Database Table)

Question : How can I manage easily a directory of mixed geometries shapefiles to deal with them depending on their geometry ? It doesn't seem as easy as with other readers.

Please feel free to add comments if I'm not clear enough.

FME version : 2015


texts:
I have a warehouse that has many storage yards. I have an order. It has items list and I have to collect them from storage yards. Storage yards are not equal each other. Their length are different from others. I have CAD drawing of the warehouse. I convert it meaningfully for spatial queries with MapInfo and post it to PostgreSQL.

I use a point to get location of storage yards and connect them each other with lines.

I need to get shortest path between warehouse enter gate and exit gate when I collect the items from storage yards.

Now here is my doubts,

Should I use which version of PostgreSQL and PostGIS?
When I research TSP I see various algorithm exists. Should I use which technique?
How can I compare these technique each other?
Should I import my points and lines to PostgreSQL as x,y coordinate as each one is one column or both are in one column?
What should I give cost for each line? Should it be 1 for every line?


texts: Does the common flow accumulation techniques still work over flat area?

texts:
My terrain is a highly developed area--which means that

It involves lots of flat area ( gradient=0) because they have been properly paved.
And it involves sudden drop, such as retaining wall, or drain drop.

These two features make me wonder whether the usual flow accumulation algorithms ( like D8) can actually work on such terrain?

I want to do streamline analysis on it, but I am quite unsure whether the usual flow accumulation algorithms ( like D8) can actually work on such terrain? Or, as my previous question frames it, is there any open source flow accumulation algorithm in C++/C# on Windows that handles these features nicely?

Note: I would have to conduct watershed analysis at modelling stage ( because we need to study the waterflow path even before the development begins, to avoid flooding at unwanted area), which means that the DEM is based on my modelling and not something out there on the actual terrain (yet).


texts: Automatically selecting feature when another feature in different layer shares attribute data?

texts:
I have a mapset displaying foraging areas, each with a unique ID. I spatially joined the mapbooks layer, so each mapbook sheet identifies the foraging area it is depicting (1:1), so sheet 5 is foraging area 18. So when sheet 5 draws, I want it to select foraging area 18, which has a matching attribute.  I want to do this so I can change the symbology of the target feature to make the map stand out more.

I'm using ArcGIS 9.3.x with the mapbooks extension.


texts: POST to the CartoDB imports API from Google AppEngine / Python

texts:
Creating tables using the CartoDB SQL api is a known bug/missing feature, CREATE TABLE queries sent via the SQL API create shadow tables that never show up in the dashboard. One potential way around this is using the new imports api to post a small csv file in the desired schema to create a table that then can be updated/inserted to via the sql api.

http://blog.cartodb.com/post/54101913823/got-files-weve-got-a-import-api

This works well with files sent via curl. It would be nice to know a little more about what parameters the imports api expects and can respond to. Specifically, I am attempting to post using urlfetch in AppEngine, and not having much luck.

Like this..

url = 'http://mycartodb.cartodb.com/api/v1/imports?api_key=myapikey'

with open('testtable.txt', 'r') as fh: #testtable is a simple 3 line csv
    payload = fh.read()

result = urlfetch.fetch(
    url=url,
    payload=payload,
    method=urlfetch.POST
).content

logging.info(result)

These requests successfully hit the api...

{
    "item_queue_id":"d18799a5-3849-4058-ac23-49ea7d6e194e",
    "success":true
}

...but when I query the job status it reports an error of unknown origin. One hitch may be that the "Content-Length" header cannot be set in AppEngine, and maybe a successful imports api POST request requires this?

{
    "id": "blah-blah-blah",
    "user_id": "blah-blah-blah",
    "table_id": null,
    "data_type": "file",
    "table_name": null,
    "state": "failure",
    "error_code": 99999,
    "queue_id": "blah-blah-blah",
    "get_error_text": {
        "title": "Unknown",
        "what_about": "Sorry, something went wrong and we're not sure what. Try uploading your file again, or contact us and we'll try to help you quickly."
    },
    "tables_created_count": null,
    "synchronization_id": null,
    "success": false
}


texts: Analyzing spatial distribution of invasive specie across landscape?

texts:
I'm studying the spatial distribution of an invasive specie introduced in a natural area.
This plant, that thrives in sunlight, grows in thickets (patches).

The population of my study site (10x10 km) was mapped using photo-interpretation.

I get polygons, each with a rate of invasion.

I'd like to understand the relative effects of several variables on its spatial distribution : open/closed area, burned areas, distance from roads and rivers, soil type and altitude.

I splatted the area into cells which have a value for each of the variables.

What kind of analysis would you do ?

I'd go for a GLM, using R, on the rate of invasion [0;X] in which the specie is present ~ variables.

Are there any other approaches ?

I work with QGIS & R.


texts: I am using QGIS 2.4 and how do I use data from a zip file.I am trying to add a layer of a base map for south africa into my workspace by am not getting any joy.


texts: Understanding terms in Length of Degree formula?

texts:
Online calculators such as http://www.csgnetwork.com/degreelenllavcalc.html (view page source) use the below formulas to get meters per degree. I understand in general how distance per degree varies depending on latitude location, but I do not understand how that translates to the below. More specifically, where do the constants, the 3 "cos" terms in each formula, and coefficients (2, 4, 6; 3, and 5) for "lat" come from?

    // Set up "Constants"
    m1 = 111132.92;     // latitude calculation term 1
    m2 = -559.82;       // latitude calculation term 2
    m3 = 1.175;         // latitude calculation term 3
    m4 = -0.0023;       // latitude calculation term 4
    p1 = 111412.84;     // longitude calculation term 1
    p2 = -93.5;         // longitude calculation term 2
    p3 = 0.118;         // longitude calculation term 3

    // Calculate the length of a degree of latitude and longitude in meters
    latlen = m1 + (m2 * Math.cos(2 * lat)) + (m3 * Math.cos(4 * lat)) +
            (m4 * Math.cos(6 * lat));
    longlen = (p1 * Math.cos(lat)) + (p2 * Math.cos(3 * lat)) +
                (p3 * Math.cos(5 * lat));


texts: I have a collection of map rasters in various projections. How do I use QGIS to convert them to simple latlong?

texts:
My name is Samuel Philippian and I am new to this forum. I am an American high school student who is interested in tectonics and geophysics. I am currently building a global database of sutures, major faults and thrusts, orogenic events, magmatic episodes and plate boundaries with the goal of building a framework of Earth's tectonic evolution over the past three billion years.

I am surfing the internet for maps to digitize and add to my database. Unfortunately, I have been mostly roadblocked for the past several months because most of the maps I am using as reference materials are in the wrong map projection. GPlates, the software I am using for this project, only accepts rasters that are in simple equirectangular, or latlong, projection. But most of the maps that I have found are in other projections in order to minimize distortion.

I emailed a geophysics researcher about this problem and showed him several of the problematic maps, asking him for help in georeferencing them. He said that he unfortunately didn't have the time to undertake that task for me, but he did recommend QGIS to me as a software that was able to georeference the images.

I downloaded QGIS3.12 and experimented with it for a while to try and accomplish a correct reprojection of a few of the images into simple latlong, but despite watching several QGIS tutorials and reading a large amount of documentation in online user manuals and forums like this one, I was never able to produce any georeferenced maps without unacceptable amounts of error.

None of these maps contain any information whatsoever about what projection they are set in. Most of them don't have a grid, only some latitude and longitude tick marks along their edges. Some of them have no latitude or longitude information at all, and show only coastlines and political boundaries for reference. I couldn't find any helpful information online about how to figure out what projection a map is in when that information is not given. I feel like I am trying to squeeze blood from a stone with this project.

I am clueless, and absolute beginner, and need help. For reference, I have attached one of the problematic images that is set in an unknown map projection. Help would be greatly appreciated!


texts: How to design a bus station query database schema?

texts:
We have bus station data, and we want to build an application which will provide the line/multiple line for a given start station and the end station.

For example user try to get a bus line suggestion from station1 to station2.

If there is a bus line which can cover both station1 and station2, this line should be returned. The result may look like this:

Step1: station1 -- station2

If there is no direct bus line between station1 and station2, then the application should try to find the exchange plan,for example, the result may look like this:

Step1: station1 -- exchangestation

Step2: exchangestation -- station2

Now we have the data, but we do not know how to design the data model in the database, how to create the schema to make the query efficient?

=============================================

Update:

For example, I have four bus lines(actually two) each with a different color:

l1: A-B-C-D
l2: D-B-A (in fact, the l2 is the reverse of l1 except the l2 will skip station `C`.
l3: E-B-D-F
l4: F-D-B-E(reverse of l3)

Now if we want to save these information in the database to query the bus line plan for given one station to another, then how many tables do we need, and what should be put into each table ?


texts: Finding location points where travel angle changes using R?

texts:
I want to use R to find out the point locations at which my vehicle (equiped with GPS) started to turn. The data contains coordinates and time information. I have sort my data based on time sequence:

 time_stamp <- as.POSIXct(paste(field$Date, field$Time), format="%d/%m/%Y %H:%M:%S")
 field <- data.frame(field,time_stamp)
 field <- field[order(field$time_stamp),]

Then I borrow the codes from post below to calculate the angles of each three sequential points:
https://stackoverflow.com/questions/11184381/calculating-angle-from-latitude-and-longitude

And modified it a little bit.

# define coordinates

coordinates(field) <- c("Longitude", "Latitude")
proj4string(field) <- CRS("+init=epsg:4326")

# calculate angles

trackAngle <- function(xy) {
angles <- abs(c(trackAzimuth(xy), 0) -
              c(0, rev(trackAzimuth(xy[nrow(xy):1, ]))))
angles <- ifelse(angles > 180, 360 - angles, angles)
angles[is.na(angles)] <- 180
angles[-c(1, length(angles))]
}

ang <- trackAngle(field@coords)

# Just to make it the same length as my dataset
# I inserted the first and the last angles to the vector

ang <- c(ang[1],ang,ang[1])

Now I am thinking to use if-else statement, says if at this point the angle is bigger than 20 degree compared the last point, then this point will be classified as the point of turning.

new_df <- data.frame(field,ang1)
new_column <- 1:nrow(new_df)

for(i in 1:nrow(new_df))
{
  angle_changes <- ang1[i+1]-ang1[i]
  if (angle_changes>20)
   {new_column[i] <- 'turn'}
  else
   {new_column[i] <- 'not_turn'}
 }
Error in if (angle_changes > 20) { :
  missing value where TRUE/FALSE needed

Is there a simpler and faster way of achieving this, without calculating all angles first?

For reproducible example, I selected the first 20 rows. The raw dataset has ~63,000 obs.

     field1
           Longitude  Latitude
     49221   175.723 -37.82931
     49222   175.723 -37.82931
     49223   175.723 -37.82931
     49230   175.723 -37.82931
     49248   175.723 -37.82931
     49272   175.723 -37.82931
     49309   175.723 -37.82932
     49324   175.723 -37.82932
     49337   175.723 -37.82932
     49347   175.723 -37.82932
     49357   175.723 -37.82932
     49367   175.723 -37.82932
     49380   175.723 -37.82932
     49389   175.723 -37.82932
     49393   175.723 -37.82932
     49400   175.723 -37.82932
     49411   175.723 -37.82932
     49417   175.723 -37.82932
     49430   175.723 -37.82932
     49434   175.723 -37.82932

I compared trip package trackAngle with the defined trackAngle. The former output 18 results NA NA 172.0851 174.9314 179.6634 172.2936 170.6603 178.0629 176.5133 178.0556 178.0443 173.4078 177.9946 179.6743 178.1043 178.5513 177.2740 177.2740,

the later output 180.0000 180.0000 172.0851 174.9314 179.6634 172.2936 170.6603 178.0629 176.5133 178.0556 178.0443 173.4078 177.9946 179.6743 178.1043 178.5513 177.2740 177.2740.