texts: QGIS ImportError: No module named 'PyQt5.QtSvg'

texts:
Recently, I've installed QGIS 2.99 on my new Ubuntu 16.04 LTS desktop:

sudo gedit /etc/apt/sources.list

deb http://qgis.org/debian-nightly xenial main

deb-src http://qgis.org/debian-nightly xenial main

sudo apt-get update

wget -O - http://qgis.org/downloads/qgis-2016.gpg.key | gpg --import

gpg --fingerprint 073D307A618E5811

gpg --export --armor 073D307A618E5811 | sudo apt-key add -

sudo apt-get update

sudo apt-get install qgis python-qgis qgis-plugin-grass

Unfortunately, the processing toolbox is missing and there is an error when launching QGIS:

Konnte Erweiterung 'processing' nicht laden

ImportError: No module named 'PyQt5.QtSvg'
Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 313, in loadPlugin
    __import__(packageName)
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/__init__.py", line 29, in
    from processing.tools.general import *              # NOQA
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/tools/general.py", line 38, in
    from processing.core.Processing import Processing
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/core/Processing.py", line 56, in
    from processing.modeler.ModelerAlgorithmProvider import ModelerAlgorithmProvider
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/modeler/ModelerAlgorithmProvider.py", line 38, in
    from processing.modeler.EditModelAction import EditModelAction
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/modeler/EditModelAction.py", line 30, in
    from processing.modeler.ModelerDialog import ModelerDialog
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/modeler/ModelerDialog.py", line 37, in
    from qgis.PyQt.QtSvg import QSvgGenerator
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/lib/python3/dist-packages/qgis/PyQt/QtSvg.py", line 1, in
    from PyQt5.QtSvg import *
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named 'PyQt5.QtSvg'

Python-Version: 3.5.2 (default, Nov 17 2016, 17:05:23) [GCC 5.4.0 20160609]
QGIS-Version: 2.99.0-Master Master, 313ec55

Python-Pfad:
/usr/share/qgis/python
/home/benuterzname/.qgis3/python
/home/benuterzname/.qgis3/python/plugins
/usr/share/qgis/python/plugins
/usr/lib/python35.zip
/usr/lib/python3.5
/usr/lib/python3.5/plat-x86_64-linux-gnu
/usr/lib/python3.5/lib-dynload
/usr/local/lib/python3.5/dist-packages
/usr/lib/python3/dist-packages
/home/benuterzname/.qgis3//python

Does anybody know how to fix this issue?


texts: GeoServer - How to add files in GeoServer VM

texts:
I would like to add files in GeoServer. GeoServer is not in local but in computer of my company.

How it possible to add files (using Image Mosaic).


texts: Connnect to PostGIS db using QGIS - when not on localhost

texts:
I"m just getting started with PostgreSQL/PostGIS db (v9.2) and I'm wanting to connect to db while not on the localhost using QGIS.  I've been able to connect successfully on localhost, however I'm not sure what I need to do to make this non localhost connection.  I'm able to ping the db server from my desktop.

In pgAdmin within the db properties I defined the following privileges:

Then in QGIS here are my connection settings (using postgres username/password, but logged in as my domain user):

I've tried searching the web for this, however all the hits are just coming up with examples of how to connect on localhost.  I'm thinking maybe I have to add my domain info as user login?

Thanks


texts: Do we get world grid with date and time of visit of Sentinel 1 and Sentinel 2?

texts:
I have an ROI and I am interested in knowing the below for it

Tile Number
Date and Time of revisit as calendar for whole year

This will enable me to specify the tile number(s) and fetch the data for those specific revisit dates. The reason for asking specific dates is , I could fill in the gaps with data from other satellite missions like landsat8 at different dates.


texts: How to add updated tiles to ArcGIS mosaic

texts:
I have several hundred aerial imagery raster tiles (in ecw format) that I have loaded into a mosaic dataset. I would now like to 'replace' or 'overlay' an older aerial imagery layer with this dataset.
Ideally I would like to merge them into one mosaic to save space and make it easier to manage.

However it is not a simple case of removing the old tiles by identifying which ones have updates (e.g SP21NE with new SP21NE) as the old aerial imagery layer has been stored as one large .ecw file and not as individual tiles in a mosaic dataset.

Could anyone suggest a work around please? Would it be best to use the updated imagery to 'clip' the old one and then merge them together in the mosaic dataset?


texts: Dynamic column header in data table from JSON array?

texts:
I want to populate data table with dynamic column header and column data,  I can populate dynamic column data successfully, but I can't achieve dynamic column...

I am using JSON array.

My html code is:

 <body>
    <table id="example" class="display" cellspacing="0" width="100%" border="1">
    </table>
 </body>

My JavaScript code is.

  <script>
  var JSONResult = '{
     "column1":"data1",
     "column2":"data2",
     "column3":"data3",
     "columnN":"dataN"}';
  var row_dtable = new Array();
  var dtable_api = $('#example').dataTable();
  $.each(JSONResult , function(key, value) {
         row_dtable.push(value);
       });
  dtable_api.api().row.add( row_dtable ).draw( false );
  </script>


texts:
Can someone explain to me what I am doing wrong here as I think I do not properly understand postgis raster bands.

I have a set of geometries that I want to merge with a raster and created a raster as a result.

I am creating the raster as follows:-

WITH RASTERS AS (
                 SELECT rast from raster_table
                ),

     GEOMETRIES AS (
                 SELECT ST_AsRaster(wkb_geometry,
                                    RASTERS.rast,
                                    {'8BUI,8BUI,8BUI'}::TEXT[],
                                    100,100,100,
                                    0,0,0
                                    ) AS raster_geom FROM geometries_table,RASTERS
                ),
     UNION AS (
       SELECT rast as rast FROM RASTERS
         UNION ALL
       SELECT raster_geom as rast FROM GEOMETRIES
     )
     SELECT ST_AsPng(ST_Union(rast)) FROM UNION

Now, this works and creates a merged raster but I am getting results as below where the color is completely out, the whole base raster becomes red.

I am sure this is to do with the Pixeltype/data/nodata values in the ST_AsRaster call, but I don't properly understand it. I am thinking that I should get these values from the base raster. This is the result I am trying to achieve, as you can see it works fine with a different base raster type:-


texts:
I would like to create a HTML file where I have a layer selector dropdown list that would be able to display one layer from a specific CartoDB visualization (this visualization has several layers in it, but I just want to display one layer at a time on the API)

I have no development skills, so I am looking for a template to do this. Anyone here knows a template that I could use to do the job?


texts: How to close polylines that mark shape of cirque glaciers?

texts:
I have mapped a lot of cirque glaciers, using polylines to mark the shape of them.

Now I need to Close these polylines, but I do not know how to do that.

I have tried converting them to polygons (which is what would be the best thing to do for me) but since my polylines are not closed, it does not work. Since I have over 70 cirques mapped, I can not close them all by adding more vertexes to each one. (That would take some time!) I am hoping there is a way to just automatically draw a straight line between each polylines two endpoints...

Can anyone help me?

I use ArcMap 10.2.1, and am a beginner.


texts: how to add arcgis 10.1 rest service in openlayers

texts: I need to add ArcGis 10.1 rest services to openlayers. Can any one share the ideas, it would be great.


texts: GeoPandas has new warning and cannot read_file()

texts:
I have written a script long time ago which worked with no problem, using GeoPandas in order to read shapefile to my jupyter notebook.
Today I have run my script and when I import geopandas I get new warning:
import geopandas as gpd

PROJ: proj_create_from_database: SQLite error on SELECT name, type,
coordinate_system_auth_name, coordinate_system_code, datum_auth_name,
datum_code, area_of_use_auth_name, area_of_use_code, text_definition,
deprecated FROM geodetic_crs WHERE auth_name = ? AND code = ?: no such
column: area_of_use_auth_name

seems like the library is being imported, yet I cannot read shapefiles.
This is how I tried to read shapefile:
shapes = gpd.read_file('shapefiles/polygons.shp')

CRSError: Invalid projection: epsg:4326: (Internal Proj Error:
proj_create: SQLite error on SELECT name, type,
coordinate_system_auth_name, coordinate_system_code, datum_auth_name,
datum_code, area_of_use_auth_name, area_of_use_code, text_definition,
deprecated FROM geodetic_crs WHERE auth_name = ? AND code = ?: no such
column: area_of_use_auth_name)

I have not changed my script and it worked for long time, so I believe something has changed with the package.  I also used to read this specific shapefile, so nothing has changed in the shapefile itself.
Has anyone faced the same issue? is something now different with GeoPandas?how can I read the shapefile?
Edit:
geopandas version : 0.9.0


texts: In PostGIS is it possible to create a View with a unique ID?

texts:
When I create a view in PostGIS is there any way to add a unique ID to that view? Just like the "gid" field in any other PostGIS table?

Edit: Sorry I should have included this in the original post. I am using PostGresql 9.0 and PostGIS 1.5.

Ando


texts: PostGIS: ST_Clip problem with area which oversteps the raster

texts:
Hello all the PostGIS wizards!

I want to programmatically clip raster with bounding box which is overstepping the border of raster. Ok, i can do it with ST_Clip, but the problem is that ST_Clip only returns the area where is raster data. I want the data from the whole bounding box area and if there is no data it should be same as nodata value (in my case it would be just zero).

Ok here is a really simple example picture i made with MS Paint just for this post:

In the st_clip i just use ST_MakePolygon for the use of clipping the raster with geometry:

ST_Clip(rast,ST_MakePolygon(ST_GeomFromText('LINESTRING(4920 70964,4920 71476,5432 71476,5432  70964,4920 70964)')),0)

So... Is there any way i can get also the 'empty area' into the clipped raster for example as a nodata value? Or is there a any way how i can get the wanted result?

I've been searching answer for this with no success so far.


texts: Seeking offline map API like Google Maps for web?

texts:
I'm looking for an offline map API for web like Google maps (because I know I can't use this offline). I currently have a program (using js and node.js) that uses the google maps api and google roads api to produce a map, and given a set of points, returns another set of points containing roads nearest to the given points, and plots this. I need to make all of this work offline.

I've been looking at offline solutions like Open Street Maps (which a lot of people suggest), but I can't seem to figure out how to plot things. I've also looked at JOSM, but is it possible to plot on JOSM using a program rather than directly on the interface?

I'm currently working on plotting roads in Taipei, so I've already downloaded MBTiles for this city, as per some suggestions I've looked at online. I'm not sure what to do with these tiles however. Is it possible to plot points/lines on these tiles, and is there a built in function (like google roads snap to roads) that finds the nearest road?

What I've tried: I have tried Open Street Maps, JOSM, Leaflet, OpenLayers, but I can't seem to be getting the hang of any of these. Which makes me think, what if these tools aren't what I'm looking for?

Question: How do I go about plotting lines on the above listed tools? What do I do with the tiles I've already downloaded? If there simply is no hope for me to do what I want to with the above listed tools, is it possible to use something else, or is my project a lost cause?


texts: How to specify a layer ID as source for a QGIS virtual layer?

texts:
When creating a virtual layer in QGIS, we can use the layer name of any existing layer in a query, like

SELECT layer1.*, layer2.*
FROM layer1
JOIN layer2 ON layer1.id = layer2.id

However layer names aren't necessarily unique. The docs say The query can use the name (or id) of existing vector layers. I assumed this would be the Identifier we can set in a layer's metadata options, but using these in a query results in a Referenced table layer_identifier in query not found! message. (Nor is it possible to use the identifier in the Source definition of an embedded layer: Query execution error on DROP TABLE IF EXISTS "table1"; CREATE VIRTUAL TABLE "table1" USING gsVLayer('spatialite','layer_identifier',UTF-8): 1 - Provider error).

We could use embedded layers, but then my problem is that we cannot use relative paths, as far as I can see (the project should be usable on different computers and OSes).

Is there a more robust mechanism than using layer names for creating a virtual layer?


texts: Can you recycle shp and shx files in shapefile fishnet (vector grid) where only dbf data changes?

texts:
I'm periodically getting data from an external partner that is bound to a very large fishnet (vector grid) in shapefile format. The fishnet itself sits at around 400~ megabytes, but the only data that ever changes is the tabular data in the database file.

So can you safely just swap the database file instead of redownloading the entire shapefile and shx-file every time you want to exchange data? Or are there circumstances where there might be changes to the index in the shapefile after saving the same geometry with new data?

It would of course be more practical to work with raster data for a fishnet this large anyway, but its sadly not an option in this case.


texts: Checking if Geoprocessing Service (async) has any jobs that are not completed yet

texts:
I am trying to find out a way to tell if a GP service has any jobs that were submitted and are not completed yet. I have noticed that one can see the status of a specific service by using the ArcGIS Server Services Directory (by using a jobid value). However, I am not interested in a specific job; my concern is only to see if there are any jobs currently started by a GP service. I need to know if a service is busy working on any job because I need to make sure that there aren't any other jobs started when running a new one. This is because the GP service writes some data which results in some locks which I want to avoid.

The current way I do that is by using Python to access the statistics page of the service: http://server:6080/arcgis/admin/services/GPServiceName.GPServer/statistics. I am using Python for server administration using the ArcGIS REST API. I have written a piece of code that uses  this page to check if there are any busy instances (which means a service is busy working on a job). If this key's value is 0 then it means the service is not working on any job and it is safe to run another one.

I consider this to be more of a workaround and looking to figure out if this information (if a GP service has any jobs running) is exposed via REST API. It seems as it doesn't, but wanted to ask others to make sure I haven't missed anything.


texts: Semi-automatic classification with GRASS modules

texts:
I'm applying between QGIS and GRASS to apply a semi-automatic classification of aerial photos.
In practice I would like to use the combination of the grass modules in the processing of QGIS; i.cluster and i.maxlink. My problem is that with i.cluster I can not create the "spectral signature" which then must be used in the second module. Why?
In practice I use two photo areas and I left (As proof) the preset parameters but at the end of the process (Although it does not give me errors), only a final report is created and no spectral signature files. Why?


texts: Applying symbology to raster layer in ArcPy without using symbology layer (layer file)?

texts:
I have a raster layer with 3 unique values (1,2,3) and I want to apply a unique color to each using RGB values or hex codes, without the use of a symbology layer (layer file).

Is it possible to create a color map or symbology layer (layer file) entirely in ArcPy?

I want my code to add the raster layer to the map with the symbology I want (as opposed to random), but I don't want to use a pre-created symbology layer (layer file) or colormap - I would like the colors to be selected entirely in ArcPy, if possible.

arcpy.MakeRasterLayer_management(cat, "cat")
arcpy.SaveToLayerFile_management("cat", "cat.lyr", "ABSOLUTE")
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
arcpy.mapping.AddLayer(df, "cat.lyr", "TOP")
# this is where I want to apply unique colors to each pixel value (pixel values are: 1,2,3)
# but I don't want to use add colormap or apply symbology from layer, unless the symbology
# layer or colormap can be created entirely in arcpy. I would like the color decisions
# to be made entirely in arcpy, without having to make a symbology layer or colormap
# in the GUI
arcpy.RefreshActiveView()
arcpy.RefreshTOC()
del mxd, df

something like this?

colormap = open("colormap.clr", "w")
colormap.write("1 70 137 102\n")
colormap.write("2 255 176 59\n")
colormap.close()


texts: Google Earth Coordinates to Google Earth Engine

texts:
I am trying to define a polygon in GEE from a set of points I have defined in Google Earth Pro.
I have converted the Google Earth given coordinates as follows: [24°27'41.42"N], [36°23'20.51"E] = [24.461506, 36.389031], but when plotting the polygon in GEE i can't seem to obtain the same shape nor centering.
The code in GEE:
    ee.Geometry.Polygon(
        [[[24.461506, 36.389031],
          [16.445217, 40.374192],
          [11.210819, 44.119544],
          [16.188453, 54.824286],
          [20.058042, 60.245078],
          [24.364033, 59.358711],
          [27.331878, 56.753617],
          [27.096642, 51.463783],
          [30.724233, 48.432361]]], null, false);

// Create a planar polygon.
var planarPolygon = ee.Geometry(polygon, null, false);

// Display the polygons by adding them to the map.
Map.centerObject(polygon);
Map.addLayer(polygon, {color: 'FF0000'}, 'geodesic polygon');
Map.addLayer(planarPolygon, {color: '000000'}, 'planar polygon');

What is my mistake?


texts: How to change events in the Oncomplete Drag function in Openlayers

texts:
below is code that works for Openlayers drag function:

drag = new OpenLayers.Control.DragFeature(polygonLayer, {
    autoActivate: true,
    onComplete: function() {alert('Are you sure this is the right postion?')}
});

Now i,m trying to change my onComplete function:

drag = new OpenLayers.Control.DragFeature(polygonLayer, {
    autoActivate: true,
    onComplete: function (event){
        var bounds = event.features[0].geometry.getBounds();
        var answer = "bottom: " + bounds.bottom  + "\n";
        answer += "left: " + bounds.left  + "\n";
        answer += "right: " + bounds.right  + "\n";
        answer += "top: " + bounds.top  + "\n";
        alert(answer);
    }
});

I must be doing something wrong as this does not work?


texts: How can we convert a .tif with .tfw .htm into a raster KML/KMZ file to open in google earth? I have a .tif file along with a .tfw and .htm files. I wanted to open it in Google Earth. Converting the .tif into a raster in KML/KMZ format will solve the problem. I was working in MATLAB mapping toolbox and did not find a solution in there. A Matlab or any opensource solution will work.


texts: How to specify projection in print composer instead of in project properties for QGIS?

texts: I am using two map objects in QGIS print composer. One of them should be in WGS 84, being an overview of North America. The other should be in MTM 5 (close-up of a region). I can specify the projection in project properties but that sets the projections for both map objects.  Is it possible to set the projection for the map object itself?  (using QGIS 2.18 but could switch to 3.0)


texts: Is there a way to get only a single output from GDALinfo?

texts:
I want to get just the Pixel size of a raster (in an automated workflow), and gdalinfo looks like the best tool for that; The trouble is that it gives a whole lot of info which I don't need.

So instead of using awk or some other logic for parsing the whole thing, is there a way of getting only the Pixel Size from gdalinfo?


texts: I have a raster image (green colour in image). I would like to find the point in this raster that is the farthest from the edge i.e. that point from which the largest number of uninterrupted cells in a given direction to the edge?


texts: ArcGIS Spatial Analyst "Sample" returns Null values

texts:
Using "sample" to join raster values onto point data, some of my points get assigned a Null raster value.

The raster does not contain Null values and definitely looks like it has a value at the location of each of these points.

Has anyone encountered this problem before?

Coincidentally, all of the points that fail to be matched to a raster cell are in a horizontal line - could this have something to do with raster tiling and if so, how would I check?


texts: How do I query for geometry inside the Arctic Circle?

texts: I have a large PostGIS database of lat/lon points.  How can I query for any points that lie above the Arctic Circle?


texts: How creat a circle in meter unit, at GCS_WGS_1984 coordinate?

texts:
My shapefile coordinate sysytem is GCS_WGS_1984.
I use the editor tool to draw a circle.
EX:Right click to enter the radius in 100 meters.
The unit seems to be degree.
I want draw a circle which unit is in meters.


texts: Geoserver: WFS/WMS query by month on Date datatype

texts:
I have a GeoServer running, using a Postgis database as storage for a workspace with a few different layers. These layers are accessible through WMS or WFS. One of the variables in this database is stored as a Date.

How should I create a ogc:filter that would allow me to select all variables for a specific month (An example could be: "select all data from June-August in 1990-1995"). I'd like this to work on both the WMS and the WFS services, but in general I think most solutions would?

I tried to use the PropertyIsLike-operator as if date was a string, but that did not work. Here is a case where I try to get all data from the month of May:

http://localhost:8090/geoserver/testWorkspace/wfs?REQUEST=GetFeature&SERVICE=WFS&VERSION=1.1.0&OUTPUTFORMAT=json&TYPENAME=elTestLayero&PROPERTYNAME=date&FILTER=
<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
    <PropertyIsLike wildCard='*' singleChar='.' escapeChar='!'>
        <PropertyName>date</PropertyName>
        <Literal>....-05-..</Literal>
    </PropertyIsLike>
</ogc:Filter>

That would return this exception:

<ows:Exception exceptionCode="NoApplicableCode">
<ows:ExceptionText>
    java.lang.RuntimeException: java.io.IOException java.io.IOException null ERROR: operator does not exist:
        date ~~ unknown Hint: No operator matches the given name and argument type(s).
        You might need to add explicit type casts. Position: 83
</ows:ExceptionText>
</ows:Exception>

It suggests that I add explicit type cast. How can I do that?

Other possible solutions I thought of:

Make a String variable which is just the same as the date variable (or I could simply add the month as a variable) in the database, and use this the way I tried to use date. This would however make the layer(s)/database slightly bigger, could it decrease the performance?
Write an more advanced sql-view for the layer with a where clause where the month(s) come from SQL view parameters. I wouldn't prefer this solution because it would require that the layer was set up in a certain way, and I'd like my software to handle different layers in a quite dynamic way.


texts: Extracting Google Earth elevation data to Excel?

texts:
I have a polyline in Google Earth.

I am able to see the elevation profile in Google Earth when I move along the polyline but is there a way to extract the lat,long and elevation data out to Excel?


texts: Use an attribute value as parameter of next tool in QGIS modeler

texts:
I am building a model where I need to use an attribute value as parameter in the next tool. I include an illustration that will hopefully make the problem clearer.

So far I have been unable to extract the attribute value outputed by the Field Calculator tool and pass it as parameter of the Buffer tool. Does anyone know how to achive this?


texts: Looking for US National Grid to Lon/Lat Converter in Python, C, or C++

texts:
I see here that there is Fortran source code for going from USNG to lon/lat: http://www.ngs.noaa.gov/PC_PROD/USNG/

Is there some python or C/C++ code that does the same?

My understanding that this is not a function provided by Proj4/pyproj/GDAL, correct?


texts: Is there a quick way to see the fields in a table via arcpy?

texts:
I frequently find myself in Python debug mode, trying to examine the fields of a layer or table.

arcpy.Describe(lyr).fields returns an array of <geoprocessing describe field object> which isn't much use.

About 10 times a day I find myself typing this:

flds = []
for fld in arcpy.Describe(lyr).fields:
    flds.append(fld.name)

Is there a faster way to obtain a readable listing of the layer or table's fields in arcpy?


texts: Positional join with multiple overlaying geometries, each containing a different attribute

texts:
I can't find any documentation on how a spatial join in QGIS is working on the following situation:
An area layer is being spatially joined with a polyline layer in order to compare attributes of both layers. But there's some cases, in which multiple lines overlap an area, each with a different attribute. My question now is: Which attribute gets selected for the spatial join?

In the screenshot you see that situation. The yellow line is the one QGIS actually took the value of and was written into the joined table. All other crossing lines have different values and seem to be ignored. How is that?


texts: Styled Layer Descriptor importing problem

texts:
I published a shapefile with GeoServer and changed that shapefile style using uDIG via WFS and then SLD file exported to some other place. My problem is while I'm trying to upload that SLD file in GeoServer I'm getting following error: "Oops something went wrong"...

Please help


texts: Algorithm for Calculating Distance Between Two Points

texts:
I'm working on an API that identifies points within a given proximity to a specific latitude and longitude.  I'm leveraging three ESRI geometry services.  The first is the buffer service that returns all locations within the proximity to the discrete lat/lon.  The second service is a query that finds all points in my feature service that are within the buffer polygon.  These two steps result in two requests to ArcGIS Server.  I'm pushing the results to a simple table that is generated from the JSON response.  I'd like to add a column to the table that shows the distance each point is from the given lat/lon.  I was able to get the data by making requests to a third geometry service - the distance service.  The problem with this is that the API might have to make a maximum of nearly 900 separate requests to the distance service.  Virtually all calls to the API will result in no more than 10 features being returned and runtime is great up to a couple hundred.  Runtime goes up to about 25 seconds for all features.  I'm really not a fan of sending 900 requests to a REST endpoint, so I did some research and came up with an equation that will determine the distance between the two points.  The equation is below.  The trouble I'm having here is that the distances are off by a relatively small degree when compared to what ESRI's distance service returns.  I've also tested using Leaflet's distance and distanceTo tools and those values are very close to what is returned from the equation below.  I'm almost certain the discrepancy is due to the Earth_Radius parameter below.  The parameter in the equation treats the radius of Earth as a constant whereas I believe ESRI has some fancy formula that adjusts the value based on lat/lon value.

So here is the question.

(1) Does anybody know where to find the equation ESRI uses to calculate distance?  When I run the API with the equation runtime is superb.

lon1 = lon1*Math.PI/180;
lat1 = lat1*Math.PI/180;
lon2 = lon2*Math.PI/180;
lat2 = lat2*Math.PI/180;
var deltaLat = lat2 - lat1;
var deltaLon = lon2 - lon1;
var a = Math.pow(Math.sin(deltaLat/2), 2) + Math.cos(lat1) * Math.cos(lat2)
     * Math.pow(Math.sin(deltaLon/2), 2);
var c = 2 * Math.asin(Math.sqrt(a));

var EARTH_RADIUS = 6371;
theDistance = c * EARTH_RADIUS * 1000;


texts: arcpy.MakeFeatureLayer_management for Feature Class within Feature Dataset with gdb

texts:
I am trying to create a layer from a feature class that exists within a feature dataset within a geodatabase. Using the interface, I can simply add data and click through the gdb and the feature data set, so I know the data does exist.

When I try to do the same in ArcPy, using this code

import arcpy
arcpy.env.workspace = r"C:\Users\mquentel\Dropbox\PLSS\GIS"
arcpy.MakeFeatureLayer_management(r"\CadNSDI\AL_CadNSDI_V2.gdb.9-1-16\AL_CadNSDI_V2.gdb\CadastralReference\PLSSFirstDivision", "Sections")

it gives out error 000732 data set does not exist or is not supported

Following Vince's advice I have run the same code without the feature data set in the path (see below) but it returns the same error.

arcpy.MakeFeatureLayer_management(r"\CadNSDI\AL_CadNSDI_V2.gdb.9-1-16\AL_CadNSDI_V2.gdb\PLSSFirstDivision", "Sections")

Using the advice from Constructing file path within Python for Feature Class within Feature Dataset?, I checked whether the file existed:

print(arcpy.Exists(r"C:/.../AL_CadNSDI_V2.gdb/CadastralReference/PLSSFirstDivision"))

where CadastralReference is the name of the feature dataset and PLSSFirstDivision is the name of the feature class as specified when checking the layer's properties and it gives out False (i.e. the data does not exist).

Here is a screenshot of the catalog as well as of the layer's properties I get if I add the data manually through the interface.

I have tried backslashes (with r), forwardslashes and also joining paths using os.join and the names of the feature dataset and class.

Where am I making the mistake?


texts: QGIS processing framework Python loading shapefile

texts:
I would like to write a Python script in the QGIS processing framework allowing the download of files from specific addresses.

Instead of starting my script with: vectorLyr = vector

I would like to have something like: vectorLyr = QgsVectorLayer("D:\...\shape_name.shp", "shape_name", "ogr")

In fact it works when I test the algorithm but fails when I want to save the program in the treatment toolbox (script considered not valid).

I can't find any solution.

Here is a simple programm showing a layer

In the first case I will have to select the layer manually:

##vectorLyr = vector
QgsMapLayerRegistry.instance().addMapLayers([vectorLyr,], True)

In this case the selection is automated:

##vectorLyr = QgsVectorLayer("D:\...\shape_name.shp", "shape_name", "ogr")
QgsMapLayerRegistry.instance().addMapLayers([vectorLyr,], True)

I have tried both options. They work fine but the 2d can't be downloaded as a treatment processing script (I have an error message).

My objective is to write a script:

which can be executed directly from the processing treatment tool box.(ex "Carte_Affectation2" in the picture)
able to download automatically a layer from a specific address. (I don't want to select it each time i run the sript)

Illustration:

Script involving a manual selection of the layer (##vectorLyr = vector):

Script involving an automatic selection of the layer (##vectorLyr = QgsVectorLayer("D:\...\shape_name.shp", "shape_name", "ogr"))


texts: Calculate minimum area where x% of raster value total is concentrated within polygons

texts:
I have a raster layer (chlorophyll concentration) and a layer of non-overlapping polygons (regions of the ocean). For each polygon, I'd like to calculate the area where the majority of chlorophyll is concentrated – e.g., the smallest area value covered by 80% of the total chlorophyll in that polygon. I cannot think of an elegant way of achieving this, however. I have access to both ArcGIS (10.1) and QGIS (2.8 Wien), so a solution with either is fine.

One crude approach is to start at a high threshold value and calculate the amount of chlorophyll above that threshold for each polygon, then turn that into a percent of the total chlorophyll for that polygon. Incrementally, I decrease the threshold until I reach 80% for each polygon and calculate the corresponding area. This is horrifically cludgey and I know that there must be a more efficient solution. Is there a better way of doing this is ArcGIS or QGIS?


texts:
I have a shapefile with projection WGS84, when I added it to GeoServer, the Native SRS is Unknown

and for Declared SRS I wrote EPSG:4326. However, the shapefile became a mess if opened in Layer Preview

Then I changed Declared SRS to EPSG:2400 (Swedish RT90), the shapefile can be shown in Layer Preview.

What could be the problem?

And if I have Google Maps as tile map, how can I overlay this shapefile on it with OpenLayers element? What's the difference between e.g. OpenLayers.xxx and ol.xxx?


texts: Examples of building-level geovisualization

texts:
I have recently stumbled upon a presentation and paper by Duncan Smith exploring geovisualization of 2001 area classification of output areas using building level resolution of urban area (London in this example).

Two questions:

Could you point me to any other examples of such applications using social/economic/health data?

What are the advantages and disadvantages of your examples?


texts: How to get NOAA web token for data access?

texts:
I've been searching the NCDC/NCEI data access site for how to get this mythical 'token' we are supposed to request to access data.  This page says nothing about it: https://www.ncei.noaa.gov/support/access-data-service-api-user-documentation

It gives examples of how to display a web page of data (no token required) but when I try to use curl I'm told half the request is missing (it's not).  I'm wondering if using a token might fix this problem, and if so, how do I get the thing?

My request and error messages are as follows:

curl https://www.ncdc.noaa.gov/access-data-service/api/v1/data?dataset=global-marine&dataTypes=WIND_DIR,WIND_SPEED&stations=AUCE&startDate=2016-01-01&endDate=2016-01-02

      testdat.txt

[1] 5991

[2] 5992

[3] 5993

[4] 5994

[2]   Done                    dataTypes=WIND_DIR,WIND_SPEED

[3]-  Done                    stations=AUCE

[4]+  Done                    startDate=2016-01-01

bvh-2:~ brianvanthull2$ {"errorMessage":"Bad Request","errorCode":400,"errors": [{"field":"startDate","message":"A

  start date is required.","value":null},
  {"field":"endDate","message":"An end date is
  required.","value":null}]}


texts: Applying the attributes of polygons to a line segment

texts:
I am trying to classify a line segment (pipeline) based off of surrounding features. I have created polygons that contain the necessary attributes (count of  total nearby houses) that overlap this line segment.

I want to apply the count of these overlapping polygons to the line segment.


texts: Is there a way to get a viz.json url from Carto Builder?

texts:
Questions for any Carto users/developers here. I am trying to find out a way to a get a viz.json url that can be used as an overlay on a Leaflet map using Carto.js API.

However, I can't find it under my Builder account dashboard anywhere. With Carto engine its still possible but with new sign-ups it seems, only Carto Builder is available.

I would like the users to be able to add the overlays themselves but if new sign-ups are only allowing Carto Builder, I don't see how they can do that.


texts:
I have several single band, 32 bit float raster datasets (fgdbr) that i want to use to build a raster stack (depth) in numpy by doing something like this dv = np.dstack((d_array, v_array))

With reference to the code snippet below,

import arcpy
d = r'P:\Sample.gdb\D_0100yr_Max'
dr  = arcpy.Raster(d)
dr.minimum # returns  0.000562695786356926
dr.maximum # returns 2.8690600395202637
dr.noDataValue # prints nothing
dr.pixelType # returns F32

d_array = arcpy.RasterToNumPyArray(dr.catalogPath, "#", "#", "#", dr.noDataValue)

d_array.max() # returns 2.86906
d_array.min() # returns -3.4028231e+38

If dr.noDataValue prints nothing, where does the RasterToNumpyArray function read the no data value from such that d_array.min() returns -3.4028231e+38 ?


texts: How can I generate and publish raster images?

texts:
I have several hundred groups of data points (in lat/lng format). I want to generate heat maps for each group using kriging or similar interpolations and then display the generated raster images on a map server. What will be the best solution? I can't use ESRI products as they're not free. I know a little bit of Python.

Thanks!


texts: Feature class, Inspection Table and Defect table - Relationship Classes

texts:
I have a file geodatabase with three elements in it.

Lighting layer (assets) consisting of point items.
Inspection table (containing inspection information for those assets)
Defects Table (linked to inspection table as parent/child and also need to be linked to the assets pulling in a range of data from the asset attributes)

My question is how to do this properly? I need the relationship classes set up so a defect can't exist without selecting its inspection record, and to create a defect the user must create it from an asset so it can pull in a variety of data automatically from the asset layer to auto populate some of the defect fields.

How would I go about setting up those relationship classes to give me the desired outcomes?

Work flow required would be:

User creates inspection record (populates data)
User clicks asset to create defect record, some data auto populates the rest he/she fills in.


texts:
I have a shapefile with city areas, but unfortunately the shapefile only contains single lines. The shapefile looks like this:

I want to convert those lines to an area-polygon shapefile in QGIS. I tried it with the "Lines to polygons" tool (Vector > Geometry Tools > Lines to polygons), but it didn't work correctly. After using this tool, the shapefile looks like this:

Is there any other way to convert the line shapefile to an area polygon shapefile in QGIS?
I am using QGIS 2.14.3 Essen.


texts:
I am able to read .las file in Python using laspy in Windows as below:

las_filename = 'abc.las'
inFile = laspy.file.File(las_filename, mode='r')

But when I put a .laz file in this code, I get the below error:

  LaspyException: Error mapping file: Laszip was not found on the system

How to directly read a .laz file as I do with .las file or if there is a way in Python to convert .laz to .las and then read it?


texts: How to add Excel-sheets with meter coordinates to QGIS?

texts: I have an Excel-sheet where the coordinates are in meters. What is the most efficient way of inputting this file to QGIS? There are thousands of files to put so please recommend some straight-forward solution.


texts: Use selected features to generate WHERE clause in ArcGIS for Portal

texts:
I have a web map in Portal that has a point layer called Work Orders.

I believe it's possible to manually select points in the layer by using the Web AppBuilder Select Widget.

Using selected points, I would like to generate a WHERE clause by listing the WONUM values in an IN CONDITION.

Example:

wonum in ('WO1754', 'WO1756', 'WO1759', 'WO1766', 'WO1770', 'WO1772', 'WO1775', 'WO1777')

From there, I would copy the WHERE clause to a separate work order management system and paste it into the system's own filtering tool.

How can I generate a WHERE clause from selected features in Portal?

(Oracle 12c)


texts: How to remove NoData areas from modis 09 data?

texts: I have mosaicked modis09 rasters that have areas of NoData values displayed in black, I have used every option possible from raster con and set null conditions in raster calculator to setting the NoData values to 0 in properties but the black spaces would not erase. Is there any other option available in ArcGIS 10.2 for solving this issue?


texts: Problem when building network with spatialite-gui

texts:
I would like to find the shortest path between two points in the following graph and corresponding table:

I am using spatialite-gui for windows version 1.2.1 and I am having a problem when building the network.
It works fine when I build it with the following parameters :

(“OK: VirtualNetwork table ‘road_net’ successfully created”)

I can then run the SQL query and get the correct result (assuming all roads are bi-directionnal).

But when I’m trying to build the network with one-way directions (parameters bellow), I get: “It’s impossible to build a Network using the given configuration; some fatal error occurred”

I don’t understand why this won’t work. ‘oneway_fromto’ and ‘oneway_tofrom’ are filled with booleans.


texts:
I have a layer of points that represent bus stops. Each point has a value (quality of service, between 0 and 100). I need to create a grid (resolution 25m) by interpolating the points values.

For each cell of the grid, I want to assign the value of the best surrounding stop (search radius of 500m), but the values of the stops should decrease linearly with the distance between them and the grid cell (-10% every 100m).

For example, a grid cell located 100m away from a stop with value 30 and 200m away from another stop with value 80 should have the value 64 (80 - 0.2*80). The bus stop with the bad value 30 has no impact on this cell.

I tried IDW interpolation but I'm not getting the result I expect. I'm using QGIS.


texts: I can't export a project (road) from AutoCAD Map3D to shapefile (QGIS). I tried also with Leica Infinity, but the only data I have exported are my benchmarks, but not the polylines of the road. I tried many times changing the Coordinates System but still nothing to do.


texts: Merging multiple tables using SQL in CartoDB

texts:
I have 5 separate tables in CartoDB each with a varying number of columns, but all containing these 4 columns (the_geom, name, type, unit_code). I want to merge them all into 1 new table and get rid of most of the columns. I figured out that UNION ALL will get rid of the columns I don't need in the end, but I can't see how to merge the tables using SQL. Whenever I try UNION ALL it doesn't want to work. I've also seen that sometimes I need to write UNION_ALL with the underscore to make it work? Here's my sample code:

SELECT the_geom, name, type, unit_code FROM points_of_interest
UNION ALL
SELECT the_geom, altname, descriptio, unitcode FROM arch_visitor_centers
UNION ALL
SELECT the_geom, name FROM arch_trailheads
UNION ALL
SELECT the_geom, bldg_name FROM arch_restrooms
UNION ALL
SELECT the_geom, name FROM arch_parking
UNION ALL
SELECT the_geom, name, FROM arch_campgrounds
UNION ALL


texts: Merge from arcpy.ListTables() producing duplicate rows?

texts:
I was trying to manually merge file geodatabase tables and was getting Null values for the output file rows and decided to write a python script to merge the files. Here's my script: import arcpy

from arcpy import env
env.workspace = "G:/US county_Climate Hazard models/NWS_County/SUMCounty_NWS_new_frq.gdb"
arcpy.env.overwriteOutput = True
tableList = arcpy.ListTables()
arcpy.Merge_management(tableList, "C:/Documents/ArcGIS/Default.gdb/Merge_sum_freq")

The code takes some time to run but runs without errors. However, when I open ArcMap back again and check for the output, I get duplicate entries in the table I made. Here's how the output looks like:

As you see every entry is duplicated. Could someone explain why is that the case and how could I fix it? Is there something going on with the way I wrote the code?

****EDIT # 1******
I also added:

print tableList

and I got:

[u'CT_09001', u'CT_09003', u'CT_09005', u'Merge_sum_freq']

(I'm not sure why I get this "u" thing)

****EDIT #2*****
The original question is now resolved thanks to Michael Miles-Stimson. I accidentally created Merge_sum_freq inside of SUMCounty_NWS_new_frq.gdb

****EDIT #3 NEW Question *****
I have a performance related question: every geodatabase table I need to merge has 1 row and 6 fields (including OBJECTID), which you could see on the screenshot I posted. When I run the above script for 10 tables, it runs in 10sec. For 30 tables - 8sec, 264 tables - ~2min. In the end I need to merge >3000 tables (every table with 1 row, so final table will have >3000 rows). Not sure if I'm extrapolating this right, but it seems like they would complete in >20min. Is this considered a good performance? If not, is there any way I could optimize anything?


texts: Is Landsat 30m resolution satellite images accurate enough to map land cover on a 40 km² watershed?

texts:
I am trying to remote sensing land cover features (bare soil, water, vegetation, cultivated areas) on a 40 km² watershed, located in the Sahel (semi-arid west africa). I am heavily interested in their spatio-temporal trends and patterns. As such, I found Landsat images to be potentially good candidates. Free available Landsat products (C1 and C2 levels) have a resolution of 30 m, i.e the surface being represented in a single pixel is close to 900 m².

I would like to know if this spatial resolution is accurate enough for my task, given the size of my study area (40 km²) ?

Is there a general rule of thumb to follow, or even better, some references discussing the topic ?

How can I account for and estimate mapping errors subsequently ?


texts:
I am wondering if there is a way, algorithm, tools, to create, from a geometry, another one that is the border with a certain thickness.

The idea is to get something similar to that :

This looks quite feasible with square or very simple shape, but with complex polygon I am a bit lost.

I am using geotools to manipulate the geometries.

I tried to scale the geometry and intersect the original and the expanded one but this does not work when the anchor point of the scale is not in the right place and I cannot figure out get it nice with all kind of shape.

Looks fine to me with geotools:

polygon.getExteriorRing().buffer(50);


texts: Ogr2ogr: ERROR 1: Unable to find driver `PostgreSQL'

texts:
I need to import geojson files into a local PostgreSQL9.6 database on OSX system.

ogr2ogr -f 'PostgreSQL' PG:"dbname=crossmap user=duccioa" "multigraph.geojson" -nln paris.multilayer

Ogr2ogr gives this error:

ERROR 1: Unable to find driver `PostgreSQL'

My GDAL2.1 library have been installed using kyngchaos binary.
I saw there are other posts on the topic but none of them resolved my problem. One issue is that I cannot figure out how to unistall the library (for example to reinstallit with HomeBrew).
The import used to work, I don't understand what went wrong at one point. The only thing I can think of is the installation of the latest version of QGIS 2.18.


texts: Export Data From Merged File to Individual

texts:
I Have a merged shape file with 'n' number of cities in it. Exporting them by selecting their names individually is a difficult thing and wasting my time.
Is their any query in QGIS to export Data from Merged file to Multiple According to their Names?


texts:
Hopefully someone can help.

I've created a TMS Layer in OpenLayers which has some local tiles these tiles are OSGB OpenData. I'm trying to get the initial zoom to be level 4 instead of 0. However OpenLayers seems to be ignoring the zoomOffset option.

Here is an example of my layer:

 var newLayer = new OpenLayers.Layer.TMS("Test Layer", "", {
        type:"png",
        getURL:overlay_getTileURL,
        alpha:true,
        setIsBaseLayer:true,
        'buffer':0,
        transitionEffect: 'resize',
        projection:new OpenLayers.Projection(projection),
        displayProjection:new OpenLayers.Projection(displayProjection),
        units:"m",
        resolutions:[400.0, 200.0, 100.0, 50.0],
        maxResolution:"auto",
        maxExtent:new OpenLayers.Bounds(0.0, 0.0, 700000.0, 1300000.0),
        zoomOffset:4
    });

  map.addLayer(newLayer);


texts: How do you find which points connect to other points

texts: I have points that I have created from vector lines. To be more specific these points are end and start points of vector lines. But what I want to do now is identify which points connect to each other without manually doing this. Is this possible?


texts:
I'm working with photos from a camera. The camera is arranged to take NDVI values instead of RGB for each pixel. I have several folders with 50 pictures in each folder as '*.asc' files the pictures.

I am using R to automatize the process, the problem is that I can't find the best way. I want to get the information of the pixels in a region, a polygon. The same for all the pictures.

I have been opening a picture as reference with ArcGIS, I draw the polygon and save it. Then I bring it to R (f(x) readshapePoly) and they I try the function 'extract('.asc*, polygon)* but I get several errors... I am new at this world also with R, and I don´t know where am I failing... I hope you can try to help me


texts:
I am using arcgis javascript api . I want to remove a recently added layer , i wrote like this

var a = map.layerIds;
var b = map.getLayer(a[2]);
map.removeLayer(b);

But I am not able to remove any layer.


texts: Map projection process for three-dimensional point cloud

texts:
I have some points given in ECEF coordinates (or LatLongAltitude based on WGS84, which is interchangable in my understanding). Currently I am working on a mapping software which allows to visualize these points in maps. The current software version orthographically projects the points onto a plane, but future versions should support a wider range of map projections like Mercator and others. (It is not yet decided which exactly but the software design should be open to extension.)

From what I have read so far, the usual process for creating maps is

project the 3D points onto the chosen reference ellipsoid which yields (2D) LatLong coordinates (in case of WGS84 ellipsoid, I could use my current LatLong and "drop" the altitude information)
transform the LatLong coordinates into map coordinates (depending on the chosen map projection)

This actually is different from what we are doing currently: By first projecting onto the reference ellipsoid, the altitude is "lost". If then the projected point is transformed (e.g. with an orthographic projection) to map coordinates, altitude has no effect. In contrast, if the 3d point is directly projected orthographically, "altitude" (wrt to any ellipsoid) has effect on the output coordinates.

I tried to illustrate this (based on the height definition in @Farid Cher's answer):

When Point P0 is directly projected (as in our case) it becomes P0p on the Projection plane, whereas by first projecting onto the ellipsoid (to P1) and then projecting on the Projection plane, it becomes P1p (which is not P0p).

My questions are these:

1.) Is there any (practically relevant) map projection which follows the approach that we have currently implemented or do they "all" adhere to the two step process described earlier? (I.e., does the Altitude have an influence on the map coordinates?)

2.) Is it safe to assume that all points with the same LatLong coordinates (wrt any practically used reference ellipsoid) are on a straight line? (I know that depending on the choice of geodetic or geocentric coordinates, these lines for different LatLong may have a common point or not, but that's not too important for me.)


texts: How to get the bounding rectangle of the rectangle drawn in the ArcGIS using ArcObjects?

texts:
I am trying to draw the rectangle in the map and get the Rectangle bounds from it. I am able to draw it but I am not getting the rectangle bounds. I have followed this url :Follow URL

And my code as :

private bool isRectangleDrawn = false;
  public void DrawRectangle2(IActiveView activeView)
  {
     if (activeView == null)
     {
        return;
  }

   IScreenDisplay screenDisplay = activeView.ScreenDisplay;
   screenDisplay.StartDrawing(screenDisplay.hDC, (System.Int16)
   esriScreenCache.esriNoScreenCache);
   IRgbColor rgbColor = new RgbColor();
   rgbColor.Red = 255;

   ESRI.ArcGIS.Display.IColor color = rgbColor;
   ISimpleFillSymbol simpleFillSymbol = new SimpleFillSymbol();
   simpleFillSymbol.Color = color;

   ISymbol symbol = simpleFillSymbol as ISymbol;
   IRubberBand rubberBand = new RubberEnvelope();
   IGeometry geometry = rubberBand.TrackNew(screenDisplay, symbol);
   screenDisplay.SetSymbol(symbol);
   screenDisplay.DrawRectangle(geometry as ESRI.ArcGIS.Geometry.IEnvelope);
  isRectangleDrawn = true;
  screenDisplay.FinishDrawing();
}

I am not sure from where I will get the Rectangle bounds of this Rectangle.


texts: "Feature type :a unknown" exception in OWSLib WFS request

texts:
I have to script as shown below and request to get info on all features included in a feature.

from owslib.wfs import WebFeatureService
wfs11 = WebFeatureService(url="https://lacunaserver.se/geoserver/ows?")
print(wfs11.identification.title)
# [operation.name for operation in wfs11.operations]

print(list(wfs11.contents))

response = wfs11.getfeature(typename='adventuremap:badplatser',
                            bbox=(10.7154412374,53.8569329668,24.2974217818,68.2833889855),
                            srsname='urn:x-ogc:def:crs:EPSG:4326')

out = open('data.html', 'wb')
out.write(bytes(response.read(), 'UTF-8'))
out.close()

However, every time I get the same error message:

"owslib.util.ServiceException: Feature type :a unknown"

What seems to be the problem in my script? I tried several WFS but they all return the same error.


texts: GeoTools errors and warnings reading .tif file

texts:
I'm trying to read my .tiff image using GeoTools. After reading I intend to perform georeferencing of the image using the four coordinates of the image. Following is my code.
public static void main(String[] args) throws IOException, FactoryException {
        AbstractGridCoverage2DReader data = new GeoTiffReader(new FileInputStream(new File("a.tif")));
        GridCoverageFactory gcf = new GridCoverageFactory();
        CoordinateReferenceSystem crs = CRS.decode("EPSG:29901");// null;//null;

        int llx = 500000;
        int lly = 105000;
        int urx = 501000;
        int ury = 106000;

        ReferencedEnvelope referencedEnvelope = new ReferencedEnvelope(llx, urx, lly, ury, crs);
        GridCoverage2D gc = gcf.create("name", (RenderedImage) data, referencedEnvelope);
        String url = "random.tif";
        File file = new File(url);
        GeoTiffWriter writer = null;
        try {
            writer = new GeoTiffWriter(file);
        } catch (IOException e) {
            e.printStackTrace();
        }
        writer.write(gc, null);
        writer.dispose();
    }

I am getting the following exceptions and warnings.
WARNING: Coordinate Reference System is not available
Exception in thread "main" org.geotools.data.DataSourceException: Raster to Model Transformation is not available for: null
    at org.geotools.gce.geotiff.GeoTiffReader.<init>(GeoTiffReader.java:279)
    at org.geotools.gce.geotiff.GeoTiffReader.<init>(GeoTiffReader.java:192)
    at IrishGridConverter.main(IrishGridConverter.java:26)
Caused by: org.geotools.data.DataSourceException: Raster to Model Transformation is not available for: null
    at org.geotools.gce.geotiff.GeoTiffReader.getHRInfo(GeoTiffReader.java:539)
    at org.geotools.gce.geotiff.GeoTiffReader.<init>(GeoTiffReader.java:265)
    ... 2 more

what could be the potential reasons, and how can I solve?


texts: QGIS: Saving a Vector Layer renames the column field title as "cen2010_1", etc

texts:
I'm trying to create a custom subset block map of the 2010 Census shapefile for Massachusetts. I've joined population data from other DBF files and am trying to save the shapefile vector layer as a new file with the added data.

When I join the files, the Attributes Table lists the data correctly with the same Column field names intact (within QGIS). When I do a "Save As..." and try to save it out as a new ESRI file, the column field names in the DBF file get renamed to "CEN2010_1" and so on. The data is fine; it's just the column names.

I didn't see many option in the Save As dialog box other than "Encoding" that might affect it. Is there something else I should be doing to keep the column field names from getting renamed?


texts: Problem with Leaflet marker collision detection

texts:
This code cannot detect marker collision. I am using a geojson file of points.

I used a javascript plugin for collision detection. Plugin source:
https://github.com/MazeMap/Leaflet.LayerGroup.Collision

My Code is given below :

<!DOCTYPE html>

    <title>Leaflet Collision LayerGroup example</title>
    <meta charset="utf-8" />
    <meta http-equiv="X-UA-Compatible" content="IE=edge" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0">

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


texts: Displaying Telstra Vector Tile Server data in QGIS 3.14.16-Pi

texts:
I am running QGIS 3.14.16-Pi and am wanting to load the layers Telstra (Australian telephone company) provides on their GIS server.
https://staging.gcm.telstra.com.au/gisserver/rest/services/Hosted/
I have tried the Vector tiles plugin without success.
I am unsure of the syntax of the URL to enter into the Vector tiles plugin.
i have tried these below, as well as removing various levels of the URL:
https://staging.gcm.telstra.com.au/gisserver/rest/services/Hosted/CATM_D_External_antenna/{z}/{x}/{y}.pbf
https://staging.gcm.telstra.com.au/gisserver/rest/services/Hosted/CATM_D_External_antenna/VectorTileServer?f=jsapi
The data displays well on a browser when clicking through the ArcGIS JavaScript link
https://staging.gcm.telstra.com.au/gisserver/rest/services/Hosted/IoT_CatM1/VectorTileServer
I am sure if it will load it is a simple format issue or use of another plugin.


texts: How to zoom to selected layers or selected features with Openlayers

texts:
I have a map with different vector layers (sources are geojson files) on it, which are grouped in on array as it follows

layers_group = [];
layers_group.push(layer_1, layer_2);

source_group = [];
source_group.push(source_1,source_2);

layers are added to the var map with the function

map.addLayer(ocean_map);
map.addLayer(layer_1);
map.addLayer(layer_2);

Layers' style is by default transparent, and it gets a different style when a specific layer or a specific feature included in the geojson file are selected.

I have hence 2 different codes for selecting a layer or a feature in the layer. My question is how to add the instruction to zoom in when the selection is completed?

My selecting codes are as follows:

Select a specific layer

 $("#region_select").change(function() {
    //it starts the function to change styles according to the selection made
      var selectedRegion = $("#region_select").val();
      var index;
      for (index = 0; index < layers_group.length; ++index){
          var name = layers_group[index].get('name');
          if(selectedRegion === name){
          layers_group[index].setVisible(true);
          document.getElementById("messageBar").innerHTML = layers_group[index].get('name');
        } else if (selectedRegion ==="all"){
            var index;
            for (index = 0; index < layers_group.length; ++index){
                layers_group[index].setVisible(true);
            }
        }
        else {
          layers_group[index].setVisible(false);
        }

      };
    });

Select a specific feature within the file

    $("#species_auto").change(function() {
    //it starts the function to change styles according to the selection made
        var selectedSpecies = $("#species_auto").val();
        console.log('selectedSpecies',selectedSpecies);
        var index;
        for (index = 0; index < source_group.length; ++index){
            source_group[index].forEachFeature(function(feat){

            //source.forEachFeature(function(feat){
                console.log(feat);
                console.log(feat.get('Species'))
                console.log(feat.get('Species').indexOf(selectedSpecies));

    //if substring(selected text) exist within fetaure property('Speices)
    //should return any value other than -1
                if (feat.get('Species').indexOf(selectedSpecies)!=-1) {
    //so set the style on each feature
                feat.setStyle(selectStyle);
                } else {
    //and if doesnt exist switch back to the deafult style
                feat.setStyle(transparentStyle);
                }
            });
        };

    });

In the first case I tried to use

var extent = layers_group[index].getSource().getExtent();
map.getView().fit(extent, map.getSize());

and in the second case

var extent = feat.getGeometry().getExtent();
map.getView().fit(extent, map.getSize());

both of them were written within the if() condition, but they didn't work. The following error was reported

TypeError: map.getView is not a function

Does somebody have a suggestion?


texts:
I would like to get a dynamic zoom linked to a drop-down menu.
For example I have a list connected to my database with all the regions and I would like that when the user choose a region that a zoom is done on the selected region.

This is the beginning of my function (in my JS) which would allow me to zoom in on a region:

// Zoom sur un pays avec une surbrillance
  function zoomToRegion(geom){

var format = new ol.format.WKT();

var region = format.readFeature(geom, {
       dataProjection: 'EPSG:4326',
       featureProjection: 'EPSG:3857'
       });

var cadre = region.getGeometry().getExtent();
var map = $('#map').data('map');
map.getView().fit(cadre, map.getSize());

And here my drop-down menu linked to my database:

<!DOCTYPE html>

         Régions
        <input list="region" type="text" placeholder=" -- Toutes -- " >
        <datalist id="region">

                    <select name="region" id="region" >
                    <?php
                    try{
                        //Tentative de connexion à la bdd
                        $bdd = new PDO('pgsql:host=localhost; dbname=TEST','postgres','');
                    }
                    catch(Exception $e){
                        // En cas d'erreur on affiche un message et on stop tout
                        die('Erreur : '.$e->getMessage());
                    }
                    $requete ='SELECT nom_reg, St_astext(geom) as geom from "GEOMETRIES"."REGION_GEOFLA_2016" order by "nom_reg" ASC';

                    $listRegion = $bdd -> query($requete);

                    foreach($listRegion as $row){
                        //echo '<option value="'.$row[0].'">'.$row[0].'</option>';
                        echo '<option value="'.$row['nom_reg'].'" onclick="zoomToRegion(\''.$row['geom'].'\')">'.$row[0].'</option>';
                    }

                    ?>
                    </select>
        </datalist>

Would you have any ideas of what is missing to achieve this feature?


texts: PyQGIS : Selected features are not Highlighted in the Map Layout

texts:
I am working on a script (plugin) that allows user to :

select a feature by attributes,

zoom to that feature,

Create a map in layout designer focused on that feature,

My problem is that the feature is not highlighted in the map layout.
I followed this topic [using QgsHighlight Class] but the result is still the same. The feature is not highlighted.
I tried the rubberband Option. Same result as QgsHighlight Class.
Note : I am using QGIS 3.18.2


texts:
How to calculate cross track distance to great cicle line and rhumb line?

For great circle line, I've tried this formula, but it is not  accurate while the distance is less than 0.1nm.

For rhumb line, I've applied the cartesian vector projection method(see: dist_Point_to_Segment). Take one point as origin and calculate vectors based on rhumb distance and bearing. Then find the closest point on rhumb line, and calculte great circle distance on it. But I'm not sure whether this method is correct for the case.


texts: WFS request works in Geoserver but server returns error

texts:
I'm trying to send a WFS request to my Geoserver but there is something wrong. Here's what I get when I select the JSONP:

http://130.239.57.16:8080/geoserver/beta/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=beta:states&maxFeatures=50&outputFormat=text%2Fjavascript

Here's the javascript:

    var owsrootUrl = 'http://130.239.57.16:8080/geoserver/ows';

var defaultParameters = {
    service : 'WFS',
    version : '1.0',
    request : 'GetFeature',
    typeName : 'beta:admin>',
    outputFormat : 'text/javascript',
    format_options : 'callback:getJson',
    SrsName : 'EPSG:4326'
};

var parameters = L.Util.extend(defaultParameters);
var URL = owsrootUrl + L.Util.getParamString(parameters);

var WFSLayer = null;
var ajax = $.ajax({
    url : URL,
    dataType : 'jsonp',
    jsonpCallback : 'getJson',
    success : function (response) {
        WFSLayer = L.geoJson(response, {
            style: function (feature) {
                return {
                    stroke: true,
                    fillColor: '#FF0000',
                    fillOpacity: 0
                };
            },
            onEachFeature: function (feature, layer) {
                popupOptions = {maxWidth: 200};
                layer.bindPopup("Popup text, access attributes with feature.properties.ATTRIBUTE_NAME"
                    ,popupOptions);
            }
        }).addTo(map);
    }
});

Nothing shows up on the web page, and Firebug shows that the request is made to the server, but the response is:

http://schemas.opengis.net/wfs/1.0.0/OGC-exception.xsd">
      Feature type beta:states> unknown

It seems like the feature type is unknown, but I don't know now to fix it.


texts: How to use real GIS tools to make skewed area maps for infographics

texts:
I'm a junior GIS tech, trying to make a map as part of an infographic. I wanted to do one where the size of polygons varied with a particular numeric attribute.

Something like this one, where each country's size is skewed by its population, but on a much smaller scale.

Does anyone know how to do this using ArcGIS?


texts: QgsRasterLayer::identify() is missing after QGis upgrade to 2.0.1

texts:
My Ubuntu has just downloaded the latest version of QGis (1.8 -> 2.01), and my plugin failed with a plenty of bugs due to implemented changes in QGis API.

I still can not figure out how do I replace the function QgsRasterLayer::identify() returning the values of the raster layer bands:

rlayer = QgsRasterLayer(rlayer_path, rlayer_basebame)
val = rlayer.identify(QgsPoint(current_point.x, current_point.y))[1]

What should I use instead of the one?


texts: Raster image won't appear in Georeferencing toolbar

texts:
I have a raster image that was georeferenced in the wrong UTM zone. When I'm in UTM 14 (correct for my location) it won't show up but when I go to UTM 15 it shows up in the layer list for the toolbar.

To remedy this I pull the image back out, delete all the extra files (.tfwx/.tix.aux.xml/.tif.over) except for the tif and pull the image in to build the pyramids in UTM 14. What I'm confused about though is it still won't show up in the Georeferencing toolbar unless I'm in UTM 15.

I figured this might be because of the mxd's caching so I even tried it again with a fresh mxd but still no luck. When the tif is the only file remaining, there is no spatial information attached to it correct? If so, then why would the mxd still be recognizing its UTM zone as UTM 15?


texts: How to manipulate mapbox symbol icons - MapboxGL.js

texts:
I have my own Mapbox style that I build in MapBox Studio. Now I want to allow my user to move poi_label symbols. How could I do that? I am following this example Create a draggable point and also I am aware of Mapbox properties. I have added my code below where I can access the layer, now how can I make this symbol draggable/movable?

     map.on('load', function () {
      var layers = map.getStyle().layers;
      layers.forEach(function (layer) {
        if (layer.type && layer.type == "symbol"
          && layer["source-layer"] && layer["source-layer"] == "poi_label") {
          layer.draggable = true;// not working
          map.on('mouseenter', layer.id, function () {
            map.getCanvas().style.cursor = 'draggable';
          });
          map.on('mouseleave', layer.id, function () {
            map.getCanvas().style.cursor = '';
          });
        }
      });

In the Create a draggable point example, the author is using manual geojson data. Is there any way to build this geojson programmatically based on the loaded part of the map?


texts: How to prepare road network for routing in QGIS or GRASS

texts:
I have been working on a road network and I want to run shortest path algorithms on it. One problem is that some lines come very close to each other, but are not connected.

Is there a way to connect these lines when they are closer to each other than some threshold?

Thanks


texts: Resizing Leaflet CircleMarkers after zoom level changes?

texts:
I am trying to find a way to have my CircleMarkers resize themselves after zooming in or out so that they scale properly with the map.  So far, I have tried firing map.InvalidateSize() on the zoomend event but that does not seem to have any effect.  Refreshing the feature layer also does not change the behavior.  You can access my map at this link (sorry for the untrusted certificate).

My ultimate goal is for the circle markers to be reduced in size as the user zooms into the map, and then increase in size as the user zooms out.

Below is my code.

    var map = L.map('map',{ zoomControl:false, attributionControl: false, center: [38.980817, -98.760076], zoom: 4});
    L.tileLayer("https://cartodb-basemaps-b.global.ssl.fastly.net/dark_all/{z}/{x}/{y}.png", {maxZoom: 20, openLegendOnLoad: false, attribution:''}).addTo(map)
    var featureLayer = new L.featureGroup().addTo(map);
    function getZoomLevel() {
        map.on('zoomend', function() {
            map.invalidateSize();
        });
        return map.getZoom();
    };
    fetch('https://craftfamily.ddns.net/data')
      .then(response => {
        return response.json()
      })
      .then(data => {
        for (var key in data.recordset) {
          //console.log(key, data.recordset[key]);
          var lat = data.recordset[key].LATITUDE;
          var lon = data.recordset[key].LONGITUDE;
          var name = data.recordset[key].NAME;
          var location = data.recordset[key].LOCATION;
          var ticker = data.recordset[key].TICKER;
          var price = data.recordset[key].PRICE;
          var change = parseFloat(data.recordset[key].CHANGE).toFixed(2);
          var marker = L.circle([lat, lon], price*10000/getZoomLevel()).bindPopup(""+name+" ("+ticker+")
$"+price + " (" + change +")
"+location+"
"+lat+", "+lon).addTo(map);
        }
      })
      .catch(err => {
      });


texts: ArcGIS API for Javascript Not Adding Layer and/or Graphic during mouse event handler?

texts:
I am trying to add a graphic to my graphics layer during a map mouse event. In the javascript console, I can see that my graphic object is created and no error is thrown when I add the graphic to my graphics layer. However, when I try to add the graphics layer to the map, nothing happens. Same thing happens when I try to do a map.graphics.add(graphic) function.  Has anyone else run into this? Or is there a workaround for this issue?

Here's my code:

app.map.on("mouse-drag",function(evt){

    if(evt.ctrlKey){  //control key must be pressed

        var polygonSymbol = new SimpleFillSymbol(); //create the symbol for the graphic
        var polygon = new Polygon(makeArc(evt));  //create the polygon the the returned geometry from makeArc()
        var graphic = new Graphic(polygon,polygonSymbol); //create the graphic objet
        console.log(graphic);
        gLayer.add(graphic); //add the graphic to the graphic layer
        app.map.addLayer(gLayer); //add the graphic layer to the map
    }

});

Also, I know that it is not being added because I can't see the html inside the map div after the event occurs, which I usually can when I add a graphic or layer.


texts: SLD filter -it is possible to write a unique filter that would filter values that contains substering “NEW”

texts:
is it possible to make a filter where I define a < Literal > as part of a < PropertyName > ?
For example:
I want do define rule for all valuse in atribute "NOTE" which contains substing "NEW"
this is part example of code for few filters:
   <Filter>
    <PropertyIsEqualTo>
      <PropertyName>NOTE</PropertyName>
      <Literal>NEW:2011</Literal>
    </PropertyIsEqualTo>
  </Filter>

   <Filter>
    <PropertyIsEqualTo>
      <PropertyName>NOTE</PropertyName>
      <Literal>NEW:2012</Literal>
    </PropertyIsEqualTo>
  </Filter>

<Filter>
    <PropertyIsEqualTo>
      <PropertyName>NOTE</PropertyName>
      <Literal>NEW:2013</Literal>
    </PropertyIsEqualTo>
  </Filter>

it is possible to write a unique filter that would filter values that contain "NEW"
I found that is workin with java regex,and I tried :
       <Filter>
        <PropertyIsEqualTo>
          <PropertyName>NOTE</PropertyName>
          <Literal>{NEW}+</Literal>
        </PropertyIsEqualTo>
      </Filter>

but its not working
additional question and a little bit more complexed:
if I have attribute values : NEW_CAR, OLD_BIKE, OLD_BUS, BUS, NEW_CAT, NEW_DOG, BUS 26,... And I want values which doesnt have substring BUS AND CAT?


texts: Check/uncheck overlays in groupedlayercontrol.js

texts:
In a Leaflet-based app I have create several overlays, grouped with leaflet.groupedlayercontrol.js, as shown in this image:

Everything works well, but I need to uncheck certain layers, allowing the user to check them if needed (the app should open with Heatmap and Geometries on, and the others off).
I have references both to the layers and to the control object. How can I check/uncheck these overlays without removing them from the map object?


texts: How to create a Polygon around a Linestring as a GeoFence in C# .NET

texts:
I'm looking for a .net library that will help to create a polygon that I can use as a Geo-fence around a route on a map.

The following image demonstrates the result I'm after. Imagine that this polygon was projected from a line string that defined a route on a map.

I've not done any GIS mapping work before, so I'm really quite lost even as where or what to look for.

UPDATE

Ok, I've had a quick look and seem to have something working. The code to generate a polygon as a buffer around the linestring is pretty simple.

var line = new LineString(points);
var bufferOp = new BufferOp(line);
var buffer = (Polygon)bufferOp.GetResultGeometry(.00001);

But what confuses me now is the distance parameter to GetResultGeomotry. In the above code I'm specifying .00001 as the distance.. this is the only way I could get a reasonable result.

What unit of measure is that?


texts: Adding button that continually adds overlays to Leaflet?

texts:
I've created this map for my upcoming road trip and I'm trying to figure out how (or if it's possible) to create a button that continually adds overlays to the leaflet map (i.e. cycling through the different days/steps in my road trip).

Example:

This is the code I am working with (each step has a variable like "d1a" (Day 1/Step A):

var route = {"color": "#FF0000", "fillOpacity": 1, "fillColor": "#FF0000", "weight": 3, "opacity": 1};

var myIcon = L.icon({iconUrl: 'circle.png', iconSize: [8, 8]})

var d1a = L.geoJSON(day01a, {onEachFeature: onEachFeature, style: route, pointToLayer: function (feature, latlng) {return L.marker(latlng, {icon: myIcon})}});
var d1a = L.geoJSON(day01b, {onEachFeature: onEachFeature, style: route, pointToLayer: function (feature, latlng) {return L.marker(latlng, {icon: myIcon})}});

var satellite   = L.tileLayer(mbUrl, {id: 'mapbox.satellite', attribution: mbAttr}),

street  = L.tileLayer(mbUrl, {id: 'mapbox.streets', attribution: mbAttr});

var map = L.map('map', {center: [37.792965, -111.485824], zoom: 5, layers: [street, d1a, d1b]});

var baseLayers = {"Satellite": satellite, "Street": street};
var overlays = {"Step 01": d1a, "Step 02": d1b};;

L.control.layers(baseLayers, overlays).addTo(map);


texts: Tool to delete all data from a shapefile

texts:
I am looking for a tool that deletes the data in the shapefile, and then copy and pastes all the features from a layer into the blank layer.

We get updates to our layer from a separate organization, but our data fields are slightly different and we have a few of our own fields. The fields that are the same have the same field name so they paste properly when you paste a layer in editor.  I would also like to eventually make this run as a timed script, so it needs to be a python friendly method. I would prefer to keep the layer the way I have it, rather than just append our own data on the end of the updated data.


texts:
I'm looking to swap from MapInfo to QGIS at my work and im trying to see if its possbile to have a centralised symbology database. eg to be able to pull our symbology from a database rather than from within the layer itself.

From what i've been able to figure out so far is that we could store the symbology for all our various layers in either an XML or a geodatabase and them import to the asset layers in question from there via postgis connection a file or url.

Has anyone had to do something like this in the past?


texts: Expanding IDW interpolation to borders of country in ArcGIS Pro?

texts:
I am new to ArcGIS.

I have a data set with air pollution concentrations for 20 locations in the Netherlands. I would like to interpolate this data with IDW, so that I can estimate the approximate concentration for every address in the Netherlands.

I used the IDW geoprocessing tool in ArcGIS Pro. My problem is, that the estimation only is made for between the data points, which leaves the parts at the borders of the country without estimation (see image).

What I would like is to expand the size of the interpolation that all of the Netherlands is captured in the estimation.

I symbolized the Netherlands in this way:

https://learn.arcgis.com/en/projects/cartographic-creations-in-arcgis-pro/lessons/add-and-symbolize-countries.htm

I tried to define the extent of the interpolation in the IDW geoprocessing tool with the shape of the country that I acquired with the method in the link (under Environments, I set Extent to: Same as layer 'The Netherlands'). But then the interpolation is spread out way too far, like 10x the surface of the Netherlands.

Does this mean I should select the country in a different way to make it work?


texts: GeoServer vector tiles TMS service how to set buffer parameter?

texts:
First, I use rest service to publish multi layers.So I want to find a place to set Default Rendering Buffer parameter or use the rest to set this buffer parameter. That I needn't set every layer one by one.

Second, I had found that use WMS service can set buffer parameter by set an explicit buffer parameter,So I want to know can set this parameter when I use TMS service? And I found that use buffer=<bufferwidth> like this can't work well: http://host/geoserver/gwc/service/tms/1.0.0/layer@EPSG:900913@pbf/{z}/{x}/{y}.pbf?buffer=100.

Even tried WMTS like this doesn't work too:http://host/geoserver/gwc/service/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&LAYER=layer&STYLE=&TILEMATRIX=EPSG:900913:{z}&TILEMATRIXSET=EPSG:900913&FORMAT=application/vnd.mapbox-vector-tile&TILECOL={x}&TILEROW={y}&VENDORPARAMS=buffer:100

So what I need is somewhere to set this buffer parameter or some methods to set this buffer parameter when I use TMS/WMTS to loading vt layer.


texts: QGIS action on layer scope - reference to layer

texts:
In QGIS I would like to create an action of type "python" on layer scope, aim of this simple action is to print out the selected feature of the layer the action is defined on:

My simple question is: How do I refer to the layer the action is defined on in the python code?


texts: Tool to start ArcGIS for Server services automatically?

texts: Is there any tools or scripts to start ArcGIS server 10.1 services automatically when the services are stop?


texts: How to find the path after "cannot find module" message when I try to use grass tools in qgis

texts: I use mac and when I try to use grass modules in qgis (brighton) I get a "cannot find module" message.


texts: Convert all KMLs within a directory to SHP, zip by name, and move to another output folder?

texts:
I have the requirement to convert a list of ~100 KML files into shapefiles and later merge and upload the polygons as one layer into ArcGIS. I'm new to GIS and Python.

I've been able to convert multiple KML files in a directory into ESRI Shapefiles (.shp, .shx, .prj, .dbf) within the same directory by installing ogr2ogr in Conda and running this in Command Line:

<strike>for %f in (*.kml) do ogr2ogr -f "ESRI shapefile" "%f"</strike>

EDIT (correct code):

for %f in (*.kml) do ogr2ogr -f "ESRI Shapefile" "%~nf.shp" "%f"

This results in the following:

How can I add onto this to move the output files to a separate folder similar to this: C:\\Users\\Desktop\\DAT Files\\Redshift Data Testing\\Output
as well as zip each group of .shp, .shx, .prj, .dbf by location (Cherry Hill, State College, etc.) similar to this:

As a follow-up question, would it be possible to create a Python program to call and perform these actions?

I will table the requirements to merge these outputs and upload into ArcGIS for now as just getting the above down would be a HUGE win for me! (I also can perform that manually for now)

The post that I used for reference was: Using ogr2ogr to convert all shapefiles in directory?


texts: Iterate through a list of polygons' centroids in QGIS

texts:
I try to calculate the azimuth between a row of polygonsheets to print the number of the last and the next sheet inside the sheet in direction of the last and the next sheet. This should help to orient inside a big Project- Atlas with many sheets.
The script (not complete):
for feature in layer.getFeatures():
    count = layer.featureCount()
    geom = feature.geometry()
    polycent = geom.centroid().asPoint()

while i < count
    azimuth = polycent[i].azimuth(polycent[i+1])

The error code is 'float' object has no attribute 'azimuth' - but my polycent[i] is no 'float' object - isn't it? How can I ask for the azimuth between the actual and the next centroid?


texts: Lines on SVG export from QGIS are rasterized

texts:
I'm working with QGIS 3.0.2-Girona on Windows 7 64bits. I have a lot of roads on a large map (30 x 40 km approximately). I am trying to export them into vectorial format to print a map. However, the lines are rasterized and I have to make a very huge zoom (1:4000 in EPSG:2154, Lambert 93). Is it because there are too many elements?

I've seen QGIS exporting huge SVG and have tried printing to pdf but obviously, it got rasterized. I will try to find a ubuntu 64 or to run it on a VM.

Is there any other suggestion? Is this a bug from Qt libs as suggested on some posts?


texts:
Does anyone have any idea about getting the map info from an ENVI header file?

This is my file.hdr as shown in the text editor:

ENVI
description = {
C:\xxxx\yyyy\Desktop\filename.}
samples = 36000
lines   = 36000
bands   = 1
header offset = 0
file type = ENVI Standard
data type = 4
interleave = bsq
byte order = 0
map info = {Geographic Lat/Lon, 1, 1, 20, 28, 0.000222222222222222, 0.000222222222222222,WGS-84}
coordinate system string = {GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]}
band names = {
Band 1}

I would like to use the "pixel easting" (20), "pixel northing" (28) and "pixel dimension" (0.000222) in a python script. Is there any Python library about that?
It seems like a simple matter, but I still can't get through it.


texts: Creating photogrametric point cloud of an object from distributed photographs, here: Making accurate 3D model of building using UAV imagery?

texts:
Is there a way of creating a photogrametric point cloud from a set of densly distributed photographs of an object?

Application:

I am trying to investigate if creating a 3D model (point cloud) of the outside of a building from UAV (Phantom 3 Pro) imagery would be a viable alternative to using a Terrestrial Laser Scanner (TLS) to create the point cloud. The UAV has a 12 megapixel camera and GPS.

The aim is to fly the UAV around the building taking images at relatively close range and creating a point cloud from dense image matching software (hopefully Pix4D). The plan is to take the usual nadir and oblique images looking down at the building but also I will try to fly at low heights and take straight on and oblique photos looking up at the windows, eaves etc to try capture as much detail as possible.


texts: How to add CSV data with no coordinates to a polygon layer in QGIS?

texts:
I would like to add attributes from an spss database (around 300 cases) to a polygon layer in QGIS. Both datasets have a column with unique IDs. I exported the spss data to a CSV file which does not contain x,y values.

The "add delimited text layer" plugin seems to require x and y values. How can I add the CSV data to the polygons?


texts: I am making contour lines for a NetCDF file for sea surface temperature and the contour lines display the first band only (bands represents different dates). Theres few 100 bands in the NetCDF, how do I make sure the contour lines match the band I want?


texts: Issue with QGIS 3.4.4 Atlas legend not showing all map items

texts:
I have a layer within my project that has been filtered to show the associated point with the associated Atlas feature:

But when in composer view and there are both point types within an Atlas preview these point types do not appear in the legend.

They are also not present in the item properties for the legend.

Do I need to be filtering the legend?


texts: What is the interface/property in ArcObjects that corresponds to a layer's dataSource property in arcpy?

texts:
I am trying to get the "Catalog Path" of a Feature Layer's data source in the table of contents of ArcMap using ArcObjects. This is very straightforward in arcpy but seemingly not so much in ArcObjects.

In arcpy, the Layer.dataSource property returns the "Catalog Path" of the layer. For example, for an SDE feature class in a feature dataset, it returns something like:

Database Connections\MyConnection.sde\SDE.MyFeatureDataset\SDE.MyFeatureClass

The help says that the dataSource property is just the workspacePath and datasetName properties combined, but this is not exactly true. For one, the workspacePath property expands the "Database Connections" shortcut, and secondly it does not include the feature dataset name. Thus workspacePath returns something like:

c:\Documents and Settings\user\Application Data\ESRI\Desktop10.0\ArcCatalog\MyConnection.sde

Incidentally this is also the value returned by IWorkspaceName.PathName.

So is there an ArcObjects property that corresponds directly to the Layer.dataSource property, or is additional logic needed to build the same "Catalog Path" from the available ArcObjects interfaces?


texts:
This question pertains to mapping arbitrary data using web map tile technology.

Most of the time the data source for a web map is spatial in nature (OSM).

I'd like to create a site map using the cartographic sense of the word.  Rather than a boring list of links I'd like to model a site map after cartography from the Metro in Washington DC.  This could be used as a sort of alternative to a standard navigation bar.

Is it possible to churn any old image/vector into a set of tiles at different zoom levels, if I were to say, draw vectors manually with inkscape?


texts: Is it possible to rotate the entire layout view, not just the data frame?

texts:
I have a project that requires the use of a CAD-like map template that is Landscape orientation, but with the title block along the right edge, with all the text in it rotated 90 degrees. Something like this:

Except that the drawing, sheet #, etc. in mine are rotated the same way the Project Name, Drawing Name, etc. are rotated in that example.

I already have this template created, but in Portrait orientation. I could rotate the data frame 90 degrees, render it to PDF, rotate the PDF page back to Lanscape, and it will look like the example above. But what I'd like to do is rotate the entire layout view so that my Portrait template is presented to me in Landscape on screen, so I can look at and manipulate the contents of the data frame in the same orientation they'll appear when printed. Rotating all the text elements and text boxes in an ArcGIS layout is a pain and usually requires lots of small adjustments to line everything up properly after it's been rotated, so I'd rather not have to create a new template and rotate everything to match the orientation above.

I don't believe ArcGIS can actually do this, but let me know if I'm wrong..


texts: How to measure Fracture in area using ArcGIS Cut and Fill?

texts: I want to measure fracture in an area.I know I have to use cut and fill tool. But I wonder how to give input for this tool? it has Input before raster surface and Input after raster surface ! what shall I enter for these 2 options?like the image below.


texts: Automatically Assign Coordinate Priority in Route Creation

texts:
I'm creating a linear reference system for a set of routes throughout a city, measuring bus performance in both directions along each corridor. I've got a set of about 50 features (although this number may increase as the analysis changes) that need routes created for each corridor. Each corridor has two feature classes assigned to it, one for each direction along the corridor. For further steps in the analysis, each corridor's direction must be a unique feature.

Rather than manually creating a route for each corridor/direction, I would like to automate the process in Modelbuilder. However, I am stuck trying to figure out a was to automate the assignment of coordinate priority. The corridors run both N-S and E-W, and am not sure how to generate the correct priority for each feature, as it only has directionality based on how the feature was drawn, rather than route measures.

I'm currently working on a Python script (using dkwiens's answer in this thread's as a guide) that will add each FC's bearing to a field, which will need to be re-classed to NW/NE/SW/SE values that match the syntax from the Create Routes tool. I have two questions:

How do I reclassify the bearing ranges to NW/NE/SW/SE values? A table join wouldn't recognize a range of values, so I'm at a loss here.

How do I then tell the Create Routes tool to pull the bearing from this value? The tool's syntax appears to require specific inputs, rather than a field value (as far as I can tell).

I am also not averse to Python, although my scripting abilities are not amazing, and the tool is going to be distributed among multiple employees, some of whom do not know any Python scripting at all.

I'm working in ArcMap Advanced 10.4.1

----------------------UPDATE------------------------------

Using the Calculate Field step and the following code, I was able to add coordinate priority values to each features table:

Code Block:

def ifBlock(BEARING):
if 0 <= BEARING <=90:
    return "UPPER_RIGHT"
elif 90.0001 <= BEARING <= 180:
    return "LOWER_RIGHT"
elif 180.0001 <= BEARING <= 270:
    return "LOWER_LEFT"
elif 270.0001 <= BEARING <= 360:
    return "UPPER_LEFT"

Expression:

ifBlock(!BEARING!)

Working in Modelbuilder, I was able to reconstruct my process and set up what looks to be the correct workflow, save for one issue. Screencapture attached:

I need to pass the output of the "Get Field Value", which has an output which matches the list of inputs for Coordinate Priority in the "Create Routes" tool. Right now the workflow ignores the "Value" output (which is set as a precondition for Create Routes) and uses the exact string input from the "Coordinate Priority" variable bubble.

I've tried using inline variables (%Value%) in the string input for the Coordinate Priority variable, but this is returned as an invalid input for the Create Routes tool.

How can I use the "Value" output from Get Field Value as the input for either my Coordinate Priority variable bubble, or as the direct input for the Create Routes tool?

I feel like I'm quite close, but am using the wrong data type, or need to feed my output into the variable using a specific tool/parameter.


texts: ArcMap 10.1 and MapProxy: what is the differense between WMS, WMTS and "Layer.OSM"

texts:
On Ubuntu I have PostgreSQL, postgis, renderd, mapnik, mapproxy.

1) I have OSM. It is rendered by Mapnik.

2) Then I have MapProxy. So I have this OSM by WMS protocol.

3) I also can have this OSM by WMTS protocol (nowdays I do not use this protocol)

On Windows I have web-program on Aptana (OpenLayers) and ArcGis Server, ArcMap with Desktop application.

So on the OpenLayers I can display this map as:

new OpenLayers.Layer.WMS("WMS", "http://192.168.36.139:8080/service?", {layers: 'osm'}); - It is WMS protocol. Here we see one picture.
new OpenLayers.Layer.OSM("OSM", "http://192.168.36.139/osm_tiles2/${z}/${x}/${y}.png", {numZoomLevels: 19}); I do not know what this protocol is. Here we see the picture as many tiles.

Also we have ArcMap and Desktop Application, upgrading by me, based on ArcObjects.

I want to display this OSM here. So I add this OSM to ArcMap by WMS protocol. GCS_WGS_1984 and 900913. What is the difference?
Setting mouse wheel zooms in ArcMap?

But ArcMap

1) Every time renders the OSM.

2) On the different scales does not display this OSM. And returns errors.

So, my questions are:

1) Do I solve this problems(with ArcMap), If I will use WMTS?

2) What is the difference between all this protocols?


texts: Obtain all elevations points from a raster DEM

texts:
I need to obtain all elevation grid points from a raster DEM.

I see the steps required as:

Read in raster DEM (Arc Info Binary Grid format).
Create subset area using Points2One (see triangle area here: http://imgur.com/aeOPBBD).
Write subset-grid elevations to a file (I can't do this yet).  I only need elevations (ie. Z only, not XY).

Any ideas for step 3?  I'm happy with getting all grid elevations or a histogram of elevations (with frequencies for (say) 100 bins).  I will later manipulate the elevation data in EXCEL or R or Fortran or Python or ???

Note: I'm very new to using GISs and am using QGIS 1.8.0 on Windows (since I found that the easiest program to use).  Thus, I'd prefer a solution using QGIS, but I'm willing to use another free GIS program if required.  Also, I'll likely have to automate/program this extraction procedure later (ie. steps 2 and 3), but would prefer a GUI-based solution first for simplicity.

Thanks, Tom.


texts: Update query adding values from multi column case statement in QGIS

texts:
I am looking to calculate a single value from four columns using a CASE WHEN statement for each of four columns. For example, column a,b,c,d have numeric values and I want to classify them and sum the total between each of the four columns.

 UPDATE table SET total = CASE
  WHEN a = 0 THEN 0
  WHEN a BETWEEN  1 AND 25 THEN 1
  WHEN a BETWEEN 26 AND 50 THEN 2
  WHEN a BETWEEN 51 AND 75 THEN 3
  ELSE 4
 END

My end goal is to do this for each column (a,b,c,d) and sum the results across the four columns in a single SQL statement.

My idea is to do a CASE WHEN for column a + CASE WHEN for column b + CASE WHEN for column c + CASE WHEN for column d. I am not sure how to structure a query for SpatiaLite to do this.

Is this the correct approach?


texts: Shapely Distance different from Geopy/Haversine

texts:
I am running a specific analysis where we use shapely to create buffers around points (store locations) and then check if other points (user locations) are present within that buffer value. When i check the distance using shapely, it turns out to be different from the distance I get from geopy. The haversine formula agrees with Geopy and a check on google maps using the measure distance function also gives around the same distance

Here is an example:

from shapely.geometry import Point, shape
from pyproj import Proj, transform
from geopy.distance import vincenty, great_circle

pt_store=Point(transform(Proj(init='EPSG:4326'),Proj(init='EPSG:3857'),-76.799614, 39.435307))

pt_user=Point(transform(Proj(init='EPSG:4326'),Proj(init='EPSG:3857'),-76.79989,39.43604))

vincenty((39.435307,-76.799614),(39.43604,-76.79989)).meters
great_circle((39.435307,-76.799614),(39.43604,-76.79989)).meters
pt_store.distance(pt_user)

Vincenty: 84.77847691521336
Great_circle: 84.90640111682812
Shapely: 110.02637304449682
Haversine formula (http://www.movable-type.co.uk/scripts/latlong.html): 84.88

Which one is right? Shapely or others?
Also, is such a big difference(~22%) expected? Or am I missing something?