texts: Splitting attribute value and inserting data into columns in QGIS using Python

texts:
I am trying to separate a string value from a field and then put it into its releated Field.

The field name are characters in the string and its value is the number.

"Field"=PT 100BC 20cd 19 is the string value

PT= 100, BC=20 cd=19.

the string is irregular.

the column names are the characters in the string and number is the value associated with it.

I have so far done this.

from qgis.core import *
import qgis.utils

inlayer=iface.activeLayer()

inlayer.selectByExpression('"TRTMTHD1"=\'PLANT\'or "TRTMTHD1"=\'NATURAL\' ')
selection=inlayer.selectedFeatures()

for s in selection:
    d=(sl.split() for sl in s)

If they do the split the result is PT, 100BC, 20CD which I do not want.

Which command do I use to split the value in field column and then update the record?

After I run this

{from qgis.core import *
import qgis.utils

layer=iface.activeLayer()
fc=layer.featureCount()
f_list=[]
for i in range(0,fc):
    feat=layer.getFeature(i)
    f_list.append(feat['SPCOMP'])

field_list=[]
for f in layer.fields():
    if len(f.name())==2:
        field_list.append(f.name())

layer.selectByExpression('"TRTMTHD1"=\'PLANT\'or "TRTMTHD1"=\'NATURAL\' ')
selection=layer.selectedFeatures()

for string in f_list:
    if string!=None:
        for fn in field_list:
            string=string.replace(fn,',' + fn)

list_of_field_values=(string.split(',') for string in f_list)
for fvl in list_of_field_values:
    if fvl !='':
        fd,vl=fvl.split(' ')
        print('field=',fd,'value=',vd)}

I get AttributeError: 'list' object has no attribute 'split'


texts: How do I calculate the distance of a segment of line that is inside a polygon?

texts:
I would like to calculate the distance of a portion of a line that is inside a polygon, as well as the portion that is outside of the polygon, based on from/to lat/long points on the line, and the lat/long points that define the polygon boundary.

Use case: A vehicle travels from one to the next. Lat/long points are recorded at each stage of the journey. Lat/long points are also provided defining the state boundary.

How would I determine the distance travelled in each state?

I am really new to this, and I am not using any GIS software. I am hoping to write a function to perform this calculation, but I would really like to understand the math and techniques used to arrive at the answer.

Any references that can be provided would also be really helpful; I don't know all the GIS terminology, so I haven't had too much success with this case on Google so far...

It has been suggested that I use a library to help perform this function, which I am totally open to.  However, I still have no idea where to begin with that the methods that should be called, or even what the name of the process is for calculating that. (I'm guessing it's not "distance of a section of a line inside and outside a polygon".  Any light that can be shed would be great!


texts: GetLegendInfo SOAP call in ArcGIS Server 10

texts:
I am generating a legend for a Silverlight map using the GetLegendInfo method in the ArcGIS 10 SOAP Web Services.

The array of MapServerLegendInfo objects that is returned does not contain annotation layers or raster layers that are in my map service.  Does anyone know how I can get those layer types to return, or if there is a good workaround?

I know that the ESRI Silverlight API v2.1 contains it's own legend control, however I am not able to use that.

Thanks,
Seth


texts:
I have created time dimension mosaic layer with *.tif files and PostgreSQL data connection.
I tried to update layer with adding new tif files to data folder and using geo.reload() from geoserver-rest python package, but no success. Hope was, GeoServer would reload data folder and add new records to database table, but nothing happened.
Is there any way to add new tiles without adding completely new store and layer


texts: gDistance to determine closest point to polygon: does changing projection change results qualitatively?

texts:
I'm using gDistance from the rgeos package for R to find the nearest point to a polygon. The projection of both points and the polygon is +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0. My code looks something like this:
min(gDistance(MyPolygon, MyPoints, byid=TRUE))

Running this gives a result and the following warnings:
Warning messages:
1: In RGEOSDistanceFunc(spgeom1, spgeom2, byid, "rgeos_distance") :
Spatial object 1 is not projected; GEOS expects planar coordinates
2: In RGEOSDistanceFunc(spgeom1, spgeom2, byid, "rgeos_distance") :
Spatial object 2 is not projected; GEOS expects planar coordinates

My understanding of this is that gDistance isn't happy with my projection (or lack thereof). The actual magnitude of the result isn't important, just that the point is indeed the closest point to the polygon. My question is will this result change if I change the projection? That is, regardless of the actual numeric value of the distance between the polygon and point, will the same point be identified as being closest if I changed the projection?

Further details
My points are data that has a global coverage over the oceans at a 1 degree resolution. My polygons are regions of the ocean typically a few degrees in length, but generally quite narrow. Nevertheless, since the points and polygons originate from different sources, the two do not necessarily line up – which is why I need to find points closest to the polygons. That said, there will always be a point within 1 degree of a polygon. My results don't need to be very accurate (i.e., ballpark estimates).


texts: Getting node way relationship from Osm2pgsql data

texts:
I would like to ask what is the best way to extract which way a node belongs to given a Postgresql databased generated by Osm2pgsql.

For example, in https://www.openstreetmap.org/node/4961890042, node_id of 4961890042 gives two ways Farmington Road (59440591) and
Farmington Road (506461412) and I would like to have similar results by performing SQL query.

I've read some materials, it is possible that a separated table can be created by defining the style file. The current style I used is default.style and it looks like query public.planet_osm_rels is not a good idea.


texts: Why do ArcGIS Area Calculations not match Official Statistics?

texts:
I am trying to calculate area of polygons using a shapefile from the world bank where countries are defined.
When however I am trying to calcuate the area, the numbers that emerge are different from the official statistics on area.
Morever, whenever I use a different projection (e.g. equal area) and recalculate the area, I always get the same results.

Can someone tell me where the mistake could come from?

Thanks!


texts: Converting EPSG:28992 to universal (Google Maps) lat/lon

texts:
I'm trying to find a way to convert EPSG:28992 to Google Maps' lat/lon using PHP. Is this possible?

I'm trying to convert XYZ: 137634.397 455715.439 0.0.

I'm using Kadaster, which is a Dutch government product. If you need any more information, please ask!


texts: Selectively install frameworks for macos QGIS?

texts:
The KyngChaos method is the one-stop shop for Mac users to install QGIS. This page has links to the critical frameworks and python modules, which he graciously packages and maintains himself. Have people have had success selectively installing some of those components but not others? Or replacing them with updates? Which packages are mission critical, and which ones can be replaced or updated?

I have had installed the KyngChaos versions of numpy (v. 1.9.2), scipy (0.16.0), and matplotlib (v. 1.4.3). I have uninstalled these versions of numpy and scipy, removed the .pth overrides in order to update via sudo pip install -U, and have had no adverse effects (yet). But QGIS does not like it if you do the same with matplotlib.


texts: Error when using proj4 in visual C++ (cannot compile project using proj4 lib)

texts:
I am using Proj4 version 4.8.0, after compiled by vc i had  proj.dll and proj.lib
But when i used proj.lib in my project vc++, i cannot compile that project.
Here is my code :

#include <iostream>
#include "proj_api.h";
using namespace std;

int main()
{
    double x = 138494.92605;
    double y = 467792.640021;

    char *srid28992 = "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.04,49.91,465.84,-1.9848,1.7439,-9.0587,4.0772 +units=m +no_defs";
    char *srid4326 = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs";

    projPJ source = pj_init_plus(srid28992);
    projPJ target = pj_init_plus(srid4326);

    if(source==NULL || target==NULL)
        return false;

    //  x *= DEG_TO_RAD;
    //  y *= DEG_TO_RAD;

    int success = pj_transform(source, target, 1, 1, &x, &y, NULL );

    x *= RAD_TO_DEG;
    y *= RAD_TO_DEG;

    cout << success << endl << x << endl << y << endl;
}

Here is message of visual studio :

  1>------ Build started: Project: Proj4ApiExam, Configuration: Debug
  Win32 ------ 1>Source.obj : error LNK2019: unresolved external symbol
  _pj_transform referenced in function _main 1>Source.obj : error LNK2019: unresolved external symbol _pj_init_plus referenced in
  function _main


texts: OpenLayers assertion error when styling the layer

texts:
I have the following code that get the search result and then lists the results that has geometry in it.
findLocation(results: any) {
    results.map((res => {
      const format = new ol_format()
      const feature = format.readFeature(res.geom, {
        dataProjection: 'EPSG:4326',
        featureProjection: 'EPSG:3857'
      })
      this.vectorSource = new ol_source_vector({
        features: [feature]
      })
      var vectorPolygons = new ol_layer_vector({
        source: this.vectorSource,
        displayInLayerSwitcher: false,
        style: this.styleMe(feature)
      });
      console.log('vectorPolygons', vectorPolygons)
      this.mapService._map.addLayer(vectorPolygons)

    })
  }
  styleMe(feature) {
    let style
    if (feature.getGeometry().getType() == "Point") {
      style = new ol_style({
        image: new ol_style_circle({
          radius: 7,
          fill: new ol_style_fill({ color: '#0099ff' }),
          stroke: new ol_style_stroke({
            color: [51, 153, 255], width: 0
          })
        })
      })
    } else {
      style = new ol_style({
        fill: new ol_style_stroke({
          color: 'transparent'
        }),
        stroke: new ol_style_stroke({
          color: [51, 153, 255],
          width: 5
        })
      })
      return style
    }
  }

I want to style every geom that has been found in result on the map but I get an Assertion error like:
AssertionError {code: 41, name: "AssertionError", message: "Assertion failed. See https://openlayers.org/en/5.1.3/doc/errors/#41 for details."}
code: 41
message: "Assertion failed. See https://openlayers.org/en/5.1.3/doc/errors/#41 for details."
name: "AssertionError"

I think its because of the style before adding the vectorLayer because when I comment the style object, everything works.


texts: Downloading polygon files, like lake or municipal borders, from OpenStreetMap in QGIS?

texts:
I have a project that requires the polygon files of a large lake in the United States.

I found Lake Michigan in QGIS with the osmSearch plugin.

I did not find a way to download this data as a vector or raster file.

Can anyone help with this?


texts: Cannot import name QtTest running unittests for QGIS plugin

texts:
I have the following imports in my unittest code:

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4 import QtCore, QtGui, QtTest

When I run my unit test from the command line I get the error:

ImportError: cannot import name QtTest

QtCore and QtGui appear to load just fine.  Can anyone help me figure out why QtTest is different?

I'm running on a Mac OS X 10.10.5.  QGIS 2.8.3.  Python 2.7.10.  I believe or at least hope is the Python bundled with QGIS.  PyQt 4.11.3.


texts: QGIS Server OGC Features API fails on PostgreSQL layers

texts:
I'm attempting to use QGIS Server 3.18.1 to serve feature data via the new "WFS3" OGC Features API. I successfully installed the server by following the manual and using the referenced training data (https://github.com/qgis/QGIS-Training-Data/archive/v2.0.zip, world.qgs project, naturalearth.sqlite dataset) I can successfully interact with data at http://localhost:8000/qgisserver/wfs3/collections?
However I have not been able to access PostgreSQL layers through the Features API and from the logs it's not clear why. The layers simply aren't returned by the /collections endpoint. The same layers are available through QGIS Server's WMS interface, so I know that they are configured correctly in the project file and that the server has access to the database.
The only clue in the log is here:
INFO Server[94]: Using configuration file path from environment: /home/qgis/projects/pg-manual.qgs
WARNING PostGIS[94]: NOTICE: row number 0 is out of range 0..-1

Is there some technical reason why QGIS Server's OGC Features API implementation is currently unable to serve layers from PostgreSQL, even though they are available through WMS? If not, how can I further debug the issue here?


texts: OpenLayers 2 mousePosition projection on the fly

texts:
I'm displaying coordinates in epsg:3857 using ol.control.MousePosition() with OSM.
I'd like to convert these on the fly in another epsg:2154.

I call in the header :

<script> src="http://cdnjs.cloudflare.com/ajax/libs/proj4js/2.2.2/proj4-src.js" type="text/javascript"></script>

Then my html :

var map = new OpenLayers.Map({
        div: "map",
        projection: "EPSG:3857",
        controls:[
                new OpenLayers.Control.LayerSwitcher({'ascending':false}),
                new OpenLayers.Control.ScaleLine({
                    'geodesic': "true"
                }),
                 new OpenLayers.Control.Navigation(),
                 new OpenLayers.Control.KeyboardDefaults(),
                 new OpenLayers.Control.MousePosition(
                     {displayProjection: new OpenLayers.Projection("EPSG:2154")}
                 )
                ]
        });

but it doesn't seem to work : coordinates are false, even negative where they shouldn't.

Any idea?


texts: Calculating numbers and doing functions on related table and display it using ArcGIS Web AppBuilder?

texts:
I have a neighborhood layer and it have a relation with a table of transactions that have been made inside that neighborhood.
I'm working on ArcGis online and web app builder.

What I want is when the user click on a neighborhood, the user will see options to chose what type of indicators he wants like Avg,max,min and he has the ability to apply any of these function on any column he wants.

For example :

Sum of all transaction_price or avg in this neighborhood

transaction_price is a column in transaction table .

I tried Related Table Charts widget but it can't do any function on any column.

Simply I want to use Charts widget but with related table ,How ? or any widget to do this ?

Is there a widget to do something like this on related tables?


texts: Debugging QuickExport Parameters in Data Interoperability Extension for Exporting to PostgreSQL?

texts:
I am using ArcGIS Desktop 10.

I tried to configure the QuickExport parameters by examining the log file after I used the QuickExport tool itself. The only things I changed from the log file are the parameters as I want the user to be able to input the parameters to connect to a postgres database after they create it and pass in a layer. Any suggestions on debugging this? I didn't get any syntax errors, but if the tool didn't work I made it throw an exception.

arcpy.QuickExport_interop(layer, POSTGIS, dbname,"RUNTIME_MACROS,""HOST,%s,PORT,%d,USER_NAME,%s,PASSWORD,%s,GENERIC_GEOMETRY,no,\
LOWERCASE_ATTRIBUTE_NAMES,Yes"",META_MACROS,""DestHOST,%s,DestPORT,%d,DestUSER_NAME,%s,DestPASSWORD,%s,DestGENERIC_GEOMETRY,no, \
DestLOWERCASE_ATTRIBUTE_NAMES,Yes"",METAFILE,POSTGIS,COORDSYS,,__FME_DATASET_IS_SOURCE_,false" %((server),(portnum),(usrName),(passWord),(server),(portnum),(usrName),(passWord)))


texts: Which formula could I use to increase the vertical exaggeration of a pointcloud?

texts:
I have a pointcloud of a very flat area! I would like to exaggerate the Z values in the same way software visualisations do - e.g ArcScene

My method would be to create a raster DEM from the Pointcloud and then use the Raster calculator to use the formula but am not sure how? Would it be some sort of exponential/logorithmic formula?


texts: Calculating ratio between rows using ArcPy module

texts:
I want to calculate the ratio between the rows (i.e., row1/row2, row2/row3 and so on..) of the 'Frequency' column based on the 'Id' column and want to assign the values to a separate column. As expected there will be no ratio for the 1st row. Here is my data-
OID Id  Frequency
0   1   50680
1   2   49740
2   3   48620
3   4   48300
4   5   48004
5   6   45600
6   7   44384
7   8   42720
8   9   41890
9   10  39700
10  11  38530
11  12  38207
12  13  38106
13  14  37303
14  15  37205
15  16  36527
16  17  35301
17  18  31674
18  19  28890
19  20  24612


texts: Is a photo of an infrared slide taken in a light box a digital Infrared image?

texts: Many years ago we received many aerial infrared images. They were developed on transparent film like a slide but much bigger - 20cm x 35cm. If I now take a digital photo of these images, are they real digital infrared images? I am trying to ascertain what level of error occurs when using conventional photography to capture analog infrared images and whether I can confidently analyse these images to compare levels of chloryphyl and landuse.


texts: adding new points with mouse click and adding them to the existing attribute table

texts:
I am quite a beginner with programming in python and I am struggling with my code. I would like to create new points(trees) on my map using mouse click and afterwards add them to the existing attribute table. I managed to locate points with mouse click however they are not added to the table, it is always one row (even I added 10 points) and the columns(ID, tree, Lat and Long) are not there as well (only ID, Shape..). I tried to connect it directly to the existing data set but then when I locate a new tree all other trees are disappearing and attribute table is containing one row as well.

     def onMouseDownMap(self, x, y, button, shift):
         #adding points to the map
         mxd = arcpy.mapping.MapDocument ("CURRENT")
         sr = mxd.activeDataFrame.spatialReference
         pntFc = r"\\abc\home\xxx\Documents\moje dokumenty\xyz\arcgis\treepoints" #new polygon
         ofc=r"\\abc\home\xxx\Documents\moje dokumenty\xyz\arcgis\trees" #existing data set
         fields=["OID@","Tree","Lat", "Long"]
         Lat="SHAPE@X"
         Long="SHAPE@Y"
         pt = arcpy.Point()
         pt.X = float(x)
         pt.Y = float(y)
         ptGeoms = [arcpy.Geometry ("POINT", pt, sr)]
         if not arcpy.Exists (pntFc):
              arcpy.CopyFeatures_management(ptGeoms, pntFc)
         else:
             with arcpy.da.InsertCursor (pntFc, fields) as iCurs: #create insert cursor
             #insert row
                for pt in ptGeoms:
                    row=iCurs.insertRow()
               arcpy.RefreshActiveView()
        print "Complete"


texts: Which approach is faster for getting all POIs from MySQL/MariaDB with PHP/Laravel

texts:
Correct me if I'm wrong.

There are three approaches to get the nearest homes, users have created in my website:

To create a table with two columns(latitude, longitude) that both of them are float and say:

$latitude = 50;
$longitude = 60;

SELECT *
FROM my_table
WHERE
(latitude <= $latitude+10 AND latitude >= $latitude-10)
AND
(longitude <= $longitude+10 AND longitude >= $longitude-10)

that 10 here means 1km for example.

In this approach we can also use harvesine formula.

To merge those columns(latitude, longitude) to one column named point as POINT type and again search each row one by one.
To categorize multiple points(the coordinates of homes users have created) as a category for one section of a country i.e. city and if a query comes with $latitude and $longitude to see the nearest homes, I will check in which category they are stored IN ORDER NOT TO search all rows but search only the section this query(coordinate) belongs to.

As I guess approach number 1 is slow because of the conditions for each row of table and again slow if I use harvesine formula.

If I use ST_Distance it seems again it's slow because again it just has lots of calculations.

But if I use approach number 3 it seems it is faster to check each section for an specific point user is than check all rows. I know how to set point for each home however I don't know how to create multiple home positions as a section maybe in another table.

BTW in new versions of MySQL and MariaDB Spatial Indexes are supported in InnoDB.

My questions:

Is approach number 1 really slow or other ST_* functions are the same as this approach to check all rows with those formulas mentioned there one by one? Which one is faster?
Does approach number 2 do something other than simple conditions to make it faster? I mean does it make any changes when using type of POINT instead of float and using ST_* functions instead of doing it myself? I want to know whether the algorithm is different.
If approach number 3 is the fastest in these three approaches, how can I categorize points in order not to search all rows in a table?
How can I use Spatial Indexes to make it as fast as possible?
If any other approaches exist and I didn't mention, could you please tell me how can I get the nearest homes just by having coordinates in MySQL/MariaDB in PHP/Laravel?


texts: Query to find the Intersection coordinates of multiple Polyon

texts:
I am using the query below to find out the intersection coordinates of 2 Polygons and it's working fine, but how can I use or modify the query to find the intersection area coordinates for multiple (or N number of) polygons?

SELECT ST_AsText(
        ST_INTERSECTION(
                       ST_GeomFromText('POLYGON((1.25 9.5,1.375 18.25,7.25 17.875,7.375 9.5,1.25 9.5))'),
                       ST_GeomFromText('POLYGON((1.25 9.5,1.375 18.25,7.25 17.875,7.375 9.5,1.25 9.5))')
                )
        );

What other alternative do I have?


texts:
Is it possible to simplify just the selected parts of a polygon? I am trying to make a simplified version of the Australian Local Government Areas (LGAs) shapefile.
I am using the v.generalise function from GRASS in QGIS, which works really well for internal boundaries. However, it often fails to simplify complex coastlines, and I'm left with hundreds of LGAs that area nicely simplified down to a few tens of points, and a handful of polygons with 10s of thousands of points, nearly all of which area along complex coastlines.
A quick solution to this problem would be if I could simply select the non-shared points along each of the complex coastlines, and use QGIS's less careful Simplify function, but it doesn't seem like it is possible to use that on only part of a polygon. I guess I could also break some features apart, simplify the relevant parts, and re-join them, but I'm not sure what would happen to my metadata then.
Anyone have any ideas?


texts: How to find the best route through multiple points?

texts:
I am using postgreSQL, postGIS and pgrouting. Now I am really confused on how am I going to find the shortest path between multiple points.

I know how to find it between two points through dijkstra but I am confused how to accomplish the same result with more than two points.

I have done some research and found this:
https://github.com/pgRouting/pgrouting/wiki/One_to_many-Dijkstra---To-review

But still I did not understand if its the right way and how to install it.


texts: ArcGIS Symbology For OSM Street Flow Directions

texts: I have downloaded OSM data in [.shp.zip] format from geofabrik.de now I need to apply Street map symbology to them. Can you please let me know how I can display the flow direction symbol in the street layer?


texts: Where can I find Open feature form under options menu?

texts:
Going through the tutorials in the documentation for QGIS and I am stuck looking for a check box in the Options menu. The check box is labeled "Open feature form" and seems to have been moved in QGIS 2.8.

There isn't anything in the identify results section


texts: Using Delete columns (Batch processing) in QGIS?

texts:
I am having problem using the Delete column (Batch processing) under Vector Table tools. I do not understand what the Field to delete column means. I have tried using any of my layer's attribute column header names (like FLAG or NOTE) but it just keeps giving the error message:

  Wrong or missing parameter value:Deleted column (row1)

My attribute table for that layer has a value of NULL for the FLAG column. I plan on using batch processing as I have several similar layers. How can I get the batch processing to delete that column?


texts:
Rendering a grid in QGIS causes problems: The outmost grid line is not drawn on different zoom-levels:

Here it is ok:

While here it is not ok:

I have tried to change these rendering-settings, but had no luck:


texts: Landsat image exported from Google Earth Engine not displaying correctly in ArcGIS

texts:
I have run a script to produce a land use classification clipped to my study site using google earth engine. However, when i export the classified image and open it in ArcGIS, the image is cut off at the bottom. Could someone kindly assist.

Below is an example of what I am talking about. The link to my script is also available below.

Script : https://code.earthengine.google.com/94fcd072ee8a9cbf06bf6142e98fe156


texts: Defining Projection For NetCDF Files Using Lat/Lon raster

texts:
I've downloaded and extracted AMSR2 sea ice concentration netCDF files from here. Then I tried using the "Define Projection" tool in ArcGIS 10.4 (on Windows 10) to define their spatial reference (and "Edit" spatial reference in ArcCatalog). Based on the documentation, I guessed it was a custom NSIDC North Pole Stereographic projection.

Unfortunately it does not overlay with my other data in that projection at all. Under "Grid Description" in the documentation, it says

  3. Grid Description

  The data is displayed in two polar stereographic maps per day, one for
  the Arctic (Arc) and one for the Antarctic (Ant). Python PYRESAMPLE is
  used to derive a 3.125km grid according to the (corner) grid coordinates
  provided by NSIDC (http://nsidc.org/data/polar_stereo/ps_grids.html).

  Polar stereographic projection with the projection plane intersecting
  the earth at 70 degree and using the Hughes ellipsoid; with the radius of
  6378.273 km and an eccentricity (e) of 0.081816153, follows a semi-minor axis
  of 6356.88944891.
  Corner coordinates northern hemisphere ('Arc'):

  lon 168.35, lat 31.00,
  lon 102.34, lat 31.38,
  lon  -9.98, lat 34.36,
  lon -80.74, lat 33.94
  Corner coordinates southern hemisphere ('Ant'):

  lon  -42.24, lat -39.25,
  lon   42.24, lat -39.25,
  lon  135.00, lat -41.46,
  lon -135.00, lat -41.46

  The values for semi-major axis, semi-minor axis, projection plane, and the >corner
  coordinates are used to calculate the pixel loacations with the PYRESAMPLE >package.

  The geographic coordinates are provided in the following netCDF-file:
  LongitudeLatitudeGrid_3.125km_Arctic.nc

The polar stereographic projection details are the same as my NSIDC projection, but I don't know what to do with these corner coordinates.

I downloaded the "LongitudeLatitudeGrid_3.125km_Arctic.nc" file. It projects neatly on top of my other NetCDFs, and I can make raster layers for latitude and longitude, but how do I use these to build the coordinate system into rasters from the NetCDFs?


texts: Select country boundary from a raster and convert this to a polygon?

texts:
I have a DEM raster file of a country. From this, I want to produce a polygon shapefile which only contains the land boundaries of the country.

My idea is to select the parts of the raster which have values greater than 0 (therefore, selecting only the land part, and excluding the LAKES and OCEANS).

The problem is how to convert this to a polygon, with all the inner lakes and the ocean excluded?

I have tried Raster>Polygonize, but I don't really want this, just the boundary of the country excluding the water bodies.


texts: Stitch together images taken at similar times from different satellites

texts:
I assume this is a very common problem, but I don't know where to start.
Given two Google Earth Engine image collections (both containing a timeseries of images and multiple scenes), how can I construct a "mosaic" with bands from both collections for matching times?
For example, to create a dataset for January 5 covering my area of interest, I may need to stitch together:

Landsat scene "x" taken at 10:00pm on Jan 5 with Sentinel scene "a" taken at 3:00am on Jan 6
Landsat scene "y" taken at 10:07pm on Jan 5 with Sentinel scene "b" taken at 3:15am on Jan 6
Landsat scene "z" taken at 10:20pm on Jan 5 with Sentinel scene "c" taken at 3:05am on Jan 6
etc.

In other words, I need to take multiple scenes from the first collection (to cover the area of interest), and then find corresponding scenes from the second collection that: (i) overlap the first set (spatially), and (ii) which were taken at similar times.
I can use a filtering method (ee.Join.saveBest) to find images that correspond in time, but this ignores the spatial aspect.
I am looking for a way to do this using Google Earth Engine.


texts: Why isn't the service layer shape properly spread out?

texts:
I created a simple service layer polygon in network analyst (ArcGIS Desktop 10.2) with no hierarchies set in the street network. Now I don't understand why is this polygon created along some particular streets?

Why isn't is properly spread out? For example, why doesn't the polygon cover the areas/streets I have marked in red? Obviously those marked areas are also within 1 km.

The impedance is 1 km and 2 km respectively.

The streets are also well connected individual line elements as you can see the junctions.


texts: Disable zoom when reached max or min zoom level

texts: Is it possible to set the map to a min/maxZoom level so that it's not possible to zoom further and the current map-view doesn't change?


texts: Adjusting the admin_level of countries in osm-bright style?

texts:
I'm using OSM Bright style

Unfortunately, this style favours more admin_level=3 countries than admin_level=2, for instance, on the map below, Transnistria is labelled instead of Moldavia.

How to change that?


texts: Is there a faster way to see if a polyline intersects itself?

texts:
I have several large polyline feature classes and I need to find how many times each polyline intersects itself (where it intersects is of no importance).  I wrote the following code, which works fine, the problem is it's been running for several hours and hasn't finished going through the first polyline.  Basically, it takes each segment along the polyline and compares it with every other segment to see if they intersect.  The first polyline in my sample data has 1488 vertices, (some have a lot more, some have fewer).  That's 1487! comparisons for just one line.  Does anyone know of a quicker way to do this?  Some of my feature classes have over 3000 polylines!

import arcpy

filepath = r'C:\projects\TriCounty.gdb\Florida_GDL_Albers'

arcpy.env.workspace = filepath

These two functions determine if the line segments intersect by looking at orientation

def ccw(a, b, c):
    return (c.Y-a.Y)*(b.X-a.X) > (b.Y-a.Y)*(c.X-a.X)

def intersect(a, b, c, d):
    return ccw(a, c, d) != ccw(b, c, d) and ccw(a, b, c) != ccw(a, b, d)

Determine spatial reference of the input file

sr = arcpy.Describe("LinesTest").spatialReference

Create a new feature class to store info

arcpy.CreateFeatureclass_management(filepath,'Bikemap_SelfInt', "POLYLINE", spatial_reference=sr)
arcpy.AddField_management('Bikemap_SelfInt', 'intCount', "DOUBLE")
arcpy.AddField_management('Bikemap_SelfInt', 'Trip_ID', "DOUBLE")
print 'Table successfully created.'

Prepare insert cursor to insert results into the newly created feature class.

Fields = ['SHAPE@', 'intCount', 'Trip_ID']
insCur = arcpy.da.InsertCursor("Bikemap_SelfInt", Fields)

Search cursor to go through the polylines of the input file.  I know there is only one set of points so I didn't have to iterate through parts.

with arcpy.da.SearchCursor("LinesTest", ["SHAPE@", "trip_id"]) as cursor:

    for row in cursor:
        segments = row[0]
        trip_id = row[1]
        intCount = 0
        lenSegs = len(row[0][0])
        print lenSegs

        for i in range(0, lenSegs - 1, 1):
            P = segments[0][i]
            Q = segments[0][i+1]
            for j in range(1, lenSegs-1, 1):
                R = segments[0][j]
                S = segments[0][j+1]

                if intersect(P, Q, R, S):
                    intCount += 1
                    print intCount
        insCur.insertRow([segments, intCount, trip_id])
del insCur

Like I said, this works great, but I don't have months (years?) to let it run!  Any ideas on how to speed up this process?


texts: OpenLayers and Mapserver-WMS-Layer disappear if the map is draged

texts:
I have a client for the WMS-Services of my mapserver. The client program is written in javascript and it uses Openlayers (v. 2.12). The map has as a base layer  "Google satellite"-Layer and the mercator projection. The WMS layers are overlayers. Some of them have as a source data - raster data and some of them have as data - vector data.
I have no problem with WMS "raster" layers. But I have a problem with WMS "vector" layers. If my Openlayers map is dragged too much to the left and the generated request to mapserver is started, the WMS "vector" layers go to be invisible (and it happens only on the lowest zoom level). I found out that the bbox in generated requests is in some way wrong because the mapserver response is empty: I tested the automatically generated request strings as links in my browser (firefox  or Chrom). I have modified the mapfile and I added the mercator projection to it. But it did not help. Here is the begin of the map file:

MAP
  NAME "africa_southafrica_admin"
  PROJECTION
    "init=epsg:4326"
    "init=epsg:3857"
    "init=epsg:900913"
  END
  STATUS ON
  TRANSPARENT ON
  UNITS METERS
  WEB
    METADATA
      "wms_srs" "epsg:4326 epsg:3857 epsg:32733 epsg:900913"
      "wms_title" "ROCEEH Southafrica"
      "wms_onlineresource" "https://www.roceeh.uni-tuebingen.de/cgi-bin/mapserv?map=mapdata/MAP_88361.map"
       "wms_enable_request" "*"
      "ms_enable_modes" "*"
      "labelcache_map_edge_buffer" "-30"
    END
  END
  MAXSIZE 10000
  SIZE  2048  2048
  EXTENT -180.0 -89.0 180.0 89.0
  FONTSET "font/fontset"
LAYER
    NAME "SRTM_90m_Screening"
    STATUS OFF
    CLASS
      STYLE
        OUTLINECOLOR 102 255 51
      END
      NAME "Satellite_SRTM_90m"
      LABEL
        TYPE BITMAP
        SIZE SMALL
        COLOR 102 255 51
      END
    END
    TOLERANCE 3
    PROJECTION
      "init=epsg:4326"
      "init=epsg:900913"
    END
    GROUP "Satellite_Data_Index"
    LABELITEM "name"
    METADATA
      "wms_title" "SRTM 90m"
      "wms_group_title" "Satellite Data Index"
      "wms_abstract" "No information available."
      "wms_include_items" "all"
      "gml_include_items" "all"
      "gml_geometries" "the_geom"
      "gml_the_geom_type" "multipolygon"
    END
  ...

Here two Requests to Mapserver: the first don't work and the second is fine:

1)
https://www.roceeh.uni-tuebingen.de/cgi-bin/mapserv?map=mapdata/MAP_88361.map&FORMAT=image%2Fpng&LAYERS=SRTM_90m_Screening&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&SRS=EPSG%3A3857&BBOX=-45940489.454273,-7333063.5634033,-18555242.460298,12831635.991646&WIDTH=1399&HEIGHT=1030

2)https://www.roceeh.uni-tuebingen.de/cgi-bin/mapserv?map=mapdata/MAP_88361.map&FORMAT=image%2Fpng&LAYERS=SRTM_90m_Screening&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&SRS=EPSG%3A3857&BBOX=-13203427.488628,-10777010.309341,14181819.505347,9387689.245708&WIDTH=1399&HEIGHT=1030

My Question is: what I can do or how can I control the automatically generated Requests to avoid the described misbehavior.

Thanks in advance, ZK.


texts:
I have an idea of something that would be nice for the QGIS Composer and I think it doesn't exists yet, despite my research. How can I suggest this feature to QGIS developpers ?

I found this that shows the request but I'm a bit stuck at this point : https://hub.qgis.org/projects/quantum-gis/issues

I suppose I'm missing something.


texts:
How many features can a user create in a single edit session (insert into versioned table)? is there a limit ?

I tried inserting 6000 features, but I started getting an error at around 3500 and the editor refused to save my edits giving the message:

  FDO error: -2147216072

Is it possible that there is a limit on how many features one can insert per session ?


texts: Creating button in ArcMap to run Python program?

texts:
I have a script that I want to run within ArcMap when a button in a toolbar is clicked, and so far I have only been able to make the script into a script tool. I need it to run as a command, not as a geoprocessing tool.

Running it as a geoprocessing tool takes much longer than when the same code is executed from the command line window. I have just started looking into ArcObjects, but I want to get started using it if it is what I'll need to use.

If anyone has any sample code, or resources for creating a button, that would be great.


texts: In Postgres how to create columns with geometry?

texts:
I have point data stored in the postgis server and I need to create a new table with their geometry values.

How can I extract their coordinates and translate them into a table. Something like:

column1 = coordinate X | column2= coordinate Y

I need this to more easily export the point data into QGIS.


texts:
Maybe (and probably) this is stupid question, but I am totaly new to FME and geometry conversion.

I need to transform WKB field from SQLServer to GML geometry using FME software.

How can I do this?


texts: Can I get value from attribute table when I click on it for using in my plugin?

texts:
I created a plugin which helps us to calculate the azimuth(bearing) between points. Here is the user interface;

When I click the points, it must fill the blanks automatically. Of course this point data has own coordinate information (X and Y). However, I cannot achieve that. Please help me!

I am using Python 3.6 - Qgis 2.99 and PyQt5


texts: How to write a multi geometries .geojson file with Shapely/Fiona?

texts:
I am very new to fiona and shapely, and trying to do something that doesn't seems to be designed : I have a list of points and a shapely.MultiPolygon (Voronoi polygons from the previous list of points), and I would like to write a single geojson file that would contain both. But from what I saw, fiona.open uses a single schema, with a defined geometry type.

So is it possible to achieve that, using fiona and shapely ?

Note : I managed to make it work using shapely and the geojson package, but my question scope is focused on shapely & fiona. If it helps anyone, here is how I did it :

features = []
features.append(geojson.Feature(geometry=my_shapely_multi_polygon, properties={}))
features.append(geojson.Feature(geometry=my_shapely_multi_point, properties={}))
feature_collection = geojson.FeatureCollection(features)

with open('myfile.geojson', 'w+', encoding='utf-8') as geojson_file :
  geojson.dump(feature_collection, geojson_file)


texts: gdal2tiles.py how to find optimal zoom level for leaflet?

texts:
I working on a project need to tackle with deal with geotiff image file.
I'm new to leaflet and geotiff image.I got an image from client and i need to overlay it to leaflet map from geotiff. I used gdal2tiles.py -p raster [image.tif] to generater tiles map. And i tried to display it on leafletmap, however I got this display

I wondering anyone have idea to figure out, how to find out the optimal zoom level for convertion from geotiff to raster files without lossing it's quality?Here is my gdalinfo from my geotiff file.

Driver: GTiff/GeoTIFF
Files: Geotiff.tif
Size is 25070, 18673
Coordinate System is:
PROJCS["WGS 84 / UTM zone 54N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",141],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32654"]]
Origin = (297884.090830000001006,3955249.192560000345111)
Pixel Size = (0.071710000000000,-0.071710000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=pix4dmapper
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  297884.091, 3955249.193) (138d45'55.76"E, 35d43'13.77"N)
Lower Left  (  297884.091, 3953910.152) (138d45'56.97"E, 35d42'30.34"N)
Upper Right (  299681.861, 3955249.193) (138d47' 7.26"E, 35d43'15.09"N)
Lower Right (  299681.861, 3953910.152) (138d47' 8.46"E, 35d42'31.66"N)
Center      (  298782.976, 3954579.672) (138d46'32.11"E, 35d42'52.71"N)
Band 1 Block=25070x1 Type=Float32, ColorInterp=Gray
  NoData Value=-10000


texts: I try building a map using CartoDB.js. The map displays a number of points. When a point is clicked, I change the layer SQL to filter some points out. What I want is that when user click on the map, outside of any point, I can handle this event to revert the map to initial state. Will it be possible?


texts: Classified Images of RandomForest-Classification look clustered

texts:
I performed a RF-classification.
I have 6 bands, ~250 trainingpoints in a shapefile, OOB estimate of error rate: 1.56%
Picture of take trainingpoints added.
My code:

setwd("D:/BA-Workspace/sw_comp_ap_aug/")

library(sp)
library(rgdal)
library(raster)
library(randomForest)

# CREATE LIST OF RASTERS
rlist=list.files(getwd(), pattern="tif$", full.names=TRUE)

# CREATE RASTER STACK
xvars <- stack(rlist)

# CREATE RASTER STACK
xvars <- stack(rlist)
x <- coordinates(xvars)[, 1]
y <- coordinates(xvars)[, 2]

x_rst <- y_rst <- xvars[[1]]
x_rst[] <- x
y_rst[] <- y

xvars <- stack(x_rst, y_rst, xvars)
names(xvars) <- c("x", "y", "fitzefa1_rot1", "fitzefa1_gr1", "fitzefa1_bl1",     "fitzefa2_rot1", "fitzefa2_gr1", "fitzefa2_bl1")
# READ Raster TRAINING DATA
sdata <- readOGR(dsn=getwd(), layer="sw_trainshape")

# ASSIGN RASTER VALUES TO TRAINING DATA
v <- as.data.frame(extract(xvars, sdata))
sdata@data = data.frame(sdata@data, v[match(rownames(sdata@data), rownames(v)),])

sdata@data  <- sdata@data[-c(5,6)]

# RUN RF MODEL
rf.mdl <- randomForest(x=sdata@data[,3:ncol(sdata@data)],     y=as.factor(sdata@data[,"Class"]),
                   ntree=501, importance=TRUE)

# CHECK ERROR CONVERGENCE
#plot(rf.mdl)
#sdata@data

# PLOT mean decrease in accuracy VARIABLE IMPORTANCE

#varImpPlot(rf.mdl, type=1)
#setOldClass(SpatialPointsDataFrame)
# PREDICT MODEL
predict(xvars, rf.mdl, filename="RfClassPred.img", type="response",
    index=1, na.rm=TRUE, progress="window", overwrite=TRUE)

Here the classified .img

The classes are obviously wrong.
Could it be due to wrong training-points?

I have 7 classes. (No excluded class for shadows or so, because I'm just trying to devide some tree-species)


texts: QGIS Layer List upgrade to V3 C++ plugin

texts:
I was upgrading a c++ plugin that I had used in the past and almost eliminated all the the usual changes according to the most current version (3.6).  The plugin used to collect current map canvas layers in a combo list box for selection.  I am getting two errors that I can't figure out what to do with. This code compiled  without problem pre-2.0.

The first happens to be about QGSMapLayer that can't be instantiated because it is abstract class  and I am getting an error
[Severity   Code    Description Project File    Line    Source  Suppression State
Error   C2259   'QgsMapLayer': cannot instantiate abstract class    tarjan_sccplugin    C:\OSGeo4W64\apps\Qt5\include\QtCore\qlist.h    214 Build   ]

This is portion of the code generating the error.

void tarjan_sccGui::populateLayers( QString const& url )
{
    cmbNetSourceLayer->clear();

QString layerName;
    QString layerType;

// fill list of layers
QMap<QString, QgsMapLayer*> mapLayers = QgsProject::instance()->mapLayers();
QMap<QString, QgsMapLayer*>::iterator layer_it = mapLayers.begin();

for (; layer_it != mapLayers.end(); ++layer_it)
{
    QgsVectorLayer* vl = dynamic_cast<QgsVectorLayer*>(layer_it.value());
    if (!vl)
        continue;
    else

     cmbNetSourceLayer->addItem(vl->name() );
}

The second error appears to be in explicit declaration for the plugin. Something has changed in how the plugins are initiated.
[Severity   Code    Description Project File    Line    Source  Suppression State
Error (active)      invalid explicit instantiation declaration  tarjan_sccplugin    d:\GIS\v3\Quantum-GIS\src\core\qgsoptionalexpression.h  82  IntelliSense    ]

No idea what must be done to eliminate this error.


texts: Well Locations in NAD27 do not overlay well locations on aerial or satellite images, i.e. Google Earth.  Computations are based on accurate locations and would like to know if it is better to place symbols for well location on Google Earth and re-compute lat/lon coordinates or use the raw data as supplied.


texts: ArcGIS Online Popups - Display Imagery date & time properties?

texts: I am displaying linked images from traffic cameras in my AGO popups. The images get updated regularly. If I download an image and look at its file info it has associated date and time properties. Is there any way to show (in the popup) the date and/or time that the images were taken?


texts: Downloading map service from ArcGIS web service as features?

texts: I imported FEMA's Advisory Base Flood Elevation (ABFE) layer. It is an ArcGIS Web Service over the internet.  I can view it on the map, but it is listed as not selectable, and I cannot export sub layers, and I cannot view the attributes table.  It seems to be locked somehow.  I need to convert this layer into a standard layer that I can clip and manipulate.


texts: Attribute error when loading rasters in QGIS TOC in plugin

texts:
I am trying to display field names dynamically inside a widget based on the number of fields in vector layer selected in combo box. When i load a raster inside QGIS TOC and a vector layer, it is raising an error as shown below

AttributeError: 'OpenlayersLayer' object has no attribute 'pendingFields'

Below is my code

    def addLayerNames(self):
     self.dlg.comboBox.clear()
     for i in reversed(range(self.dlg.verticalLayout.count())):
        self.dlg.verticalLayout.itemAt(i).widget().setParent(None)
     layers = QgsMapLayerRegistry.instance().mapLayers().values()
     layer_list = []
      for layer in layers:
        if layer.type() == QgsMapLayer.VectorLayer:
            layer_list.append(layer.name())
      self.dlg.comboBox.addItems(layer_list)
      selectedLayerIndex = self.dlg.comboBox.currentIndex()
      print selectedLayerIndex
      selectedLayer = layers[selectedLayerIndex]
      print selectedLayer
      fields = selectedLayer.pendingFields()
      fieldnames = [field.name() for field in fields]
      for i,field in enumerate(fieldnames):
        self.labels = QLabel()
        self.linedit = QLineEdit()
        self.labels.setText(field)
        #self.flay = QFormLayout()
        #self.le = QLineEdit()
        self.linedit.setFixedWidth(200)
        self.dlg.verticalLayout.addRow(field, self.linedit)
        #self.linedit.setFixedWidth(100)
        #self.dlg.verticalLayout.addLayout(self.flay)
        #self.dlg.verticalLayout.addWidget(self.linedit)
        request = QgsFeatureRequest().setFilterFid(0)
        layer.setSelectedFeatures([0])
        self.iface.actionZoomToSelected().trigger()
        feat = layer.getFeatures(request).next()
        self.linedit.setText(str(feat[i]))
        self.linedit.setObjectName(field)
        self.dlg.pushButton.setEnabled(False)
        self.lineedits.append(self.linedit)


texts:
I am trying to rotate a viewport object in R. However, I end up getting boundaries for the rotated object.

Any idea as to how to rotate but not end up with these boxes?

Below is the code. m_cont, m_AK and m_HI are the maps.

library("grid")

vp_AK <- viewport(x = 0.15, y = 0.25, width = 0.3, height = 0.3)
vp_HI <- viewport(x = 0.36, y = 0.2, width = 0.2, height = 0.13)

tmap_mode("plot")
m_cont
print(m_AK, vp = vp_AK)
print(m_HI, vp = vp_HI)


texts:
I'm looking for a way to convert a GeoJSON File (from overpass-turbo.eu) directly into a Feature Class in an Esri geodatabase.

Normally I'm doing this by loading the GeoJSON into QGIS, Export it to a shapefile and then import it into ArcGIS Pro. There I can just export the data into my geodatabase.  The problem is the limited space within the Attribute Table of shapefiles, which leads to a loss of data which I need in my current project.

Unfortunately I don't have access to the ArcGIS Data Interoperability Tool, so I have to find an open source way.


texts:
So here is a script that is used to extract a specific file from multiple folders and subfolders. It run smooth, until the end.

import os, fnmatch, arcpy

arcpy.env.workspace = r'C:\Users\sysadmin\Desktop\DataMagement\Comp_3\SA_MGCP3.gdb'

fcMerged = 'MergedFeatureClass'
fcClipped = 'Extraction_Mine'
fcToClip = r'C:\Users\sysadmin\Desktop\DataMagement\SA.shp

filesMatched = list()

for root, subfolders, files in os.walk (r'C:\Users\sysadmin\Desktop\DataMagement\MGCP_Raw'):

    for f in fnmatch.filer(files, 'AAA010.shp'):

        filesMatched.append(os.path.join(root,f))

Here is the error

    arcpy.Merge_managment(filesMatched, fcMerged)

    arcpy.Clip_analysis(fcMerged, fcToClip, Extraction_Mine)

The error code reads:

  Traceback (most recent call last)' File
  "C:\User\sysadmin\Desktop\test.py", line 17 in < module 7 >
  arcpy.Merge_management (filesMatched, fcMerged) Attribute error:
  'Module' object has no attribute 'Merge_Management'


texts:
I have tracks (just waypoints) from my GPS in gpx format. I also have orthophoto pictures in tif format, that I would like to make a layer with
I use QGIS (1.8.0)
actually i know how to put it in, but it shows wrong
(coordinate system = WGS 84/EPSG 4326 - europe)

my goal is to put together all waypoints and export as a shp file.


texts: Changing the default symbology field when adding DWG to ArcMap 9.3.1

texts: When I add a dwg to ArcMap 9.3.1, the symbology is set, by default, to something the program calls "CAD unique entity values"(CAD Renderer?) - a combination of "Linetype", "Color" and "LineWt" fields. Is there any way to change the default so that any random DWG I load will be symbolized according to the "Layer" field (for example) and not those three field mentioned above?


texts: I am working in ArcMap 10.2.2.  I am trying to project a NED 3m dataset from a Geographic Coordinate System to a Projected Coordinate System.  The output rasters pixels do not overlap the original rasters.  I have used the snap to raster in Environments, both with the Background processing enabled and not.  I have tried running the projection in ArcCatolog and Arcmap.  Does anyone have any idea how to align these rasters.


texts: Distance between two points in a line feature

texts:
I have two layers, one with lines and the other with points. Layer one has all the points (A to J) as shown below and layer two has all the lines connecting these points. Is there an automated way (in ArcGIS or shapely or anything) to calculate the distance between every pair of points where my definition of distance is the actual length of the connecting lines?

Please note that I am not looking for line-of-sight distance between points which can be computed using Haversine or Vincenty.


texts: Exporting DWG from ArcMap need to unify segments to polylines

texts:
I am attempting to export to CAD a contour shapefile I have derived from LIDAR imagery.

Once exported I bring the contour DWG into ArchiCAD, drag ' drop, boom it opens, except:

Any curves are broken down into thousands of single small co-joined lines, and unlike Autocad, I don't have an effective-efficient way to polyline edit, make continuous sections into single polylines or splines as they'd be referred to in ArchiCAD.

Is there a way to set up my shapefile, edit or geoprocess it prior to export?


texts: Shapefiles - polygon type - is it in fact multipolygon?

texts:
I'm wondering about the ESRI shapefile with shape type 5.

When reading the spec I see that "A polygon may contain multiple
outer rings.". Does this mean that these outer rings must be combined to form a polygon or does it mean that it is in fact a "MultiPolygon" type? So that multiple polygons can be in one entry?

EDIT: Background => I'm writing a reader for shapefiles.


texts: Embedding YouTube videos in CartoDB infowindows?

texts:
I need a tutorial for how to embed YouTube videos in CartoDB Editor infowindows.

Most help I can find about this is about iframes which were apparently disabled for security reasons.

What do I have to do, step by step, to make an infowindow display a YouTube video in CartoDB Editor?

What is the HTML code I have to put in the infowindow section of the CartoDB web interface?

Is this possible without hosting the map on my own server? How do I access the .js script if I am using the CartoDB Editor web interface?

This, from the aforementioned link, does not work:

  x


texts: Cropping elements from one layer to an element from another

texts:
I'm trying to make a hex map of a country (Russia, in this case). To do this, I've opened a map of Russia in QGIS (from Natural Earth), set up my projection, and created a hex grid with MMQGIS. So far, so good.

Now I want to crop the hexes with the shape of Russia. I've tried using the intersect option, but I get this error:

"Input CRS error: One or more input layers missing coordinate
 reference information, results may not be as expected."

and an empty shapefile for my troubles. What am I doing wrong?

UPDATE: Ok, so now I have them both in the same coordinate system (based on an Albers projection using meters), but I'm still not able to get an intersection to work:

  Warnings:
  Input CRS error: One or more input layers missing coordinate reference information, results may not be as expected.
  GEOS geoprocessing error: One or more input features have invalid geometry.
  Some output geometries may be missing or invalid.

The output shapefile is empty.

I was able to generate a version that almost works using WGS84 coordinates, but that had some problems with the 180° meridian...


texts:
I have a sql script that creates several tables, a trigger, and a trigger function.  All the commands work in pgAdmin when I execute them in the "Query" section.  My goal is to put all these commands in a script so that I can set up the database on any machine relatively easily and quickly.  When I attempt to run the script through psql, however, I get the following errors on all of the lines where I attempt to use Geography types:

ERROR:  type "geography" does not exist

It doesn't make sense because I can copy and paste the exact same queries into the query editor of pgAdmin and it works just fine.  Does psql not support PostGIS?  Or maybe just not the most recent version of PostGIS, 1.5.3, which includes the new Geography types?  Is there something I can do to fix this?

An example of one of my queries:

CREATE TABLE source_imagery (
    id SERIAL PRIMARY KEY,
    image_type VARCHAR(1000),
    image_path VARCHAR(1000),
    boundary GEOGRAPHY(POLYGON, 4326),
    image_time TIMESTAMP,
    catalog_time TIMESTAMP DEFAULT CURRENT_TIMESTAMP
);

Error message:

ERROR:  type "geography" does not exist
LINE 5:     boundary GEOGRAPHY(POLYGON, 4326),
                     ^


texts: Any way to acces GDAL console in processign modeler in QGIS?

texts:
Is there a possibilty to use the GDAL console in a processing model?

As first step I'd like to resize cellsize of a raster and did it (without modeler) using the -tr X Ycommand.

Maybe there are other options which I don't know.


texts:
I need advice on how to complete the following task:

I have aeronautical charts covering parts of the US
(example: http://aeronav.faa.gov/content/aeronav/sectional_files/Albuquerque_89.zip - Each zip contains ~2 files covering both the northern and the southern part of the area).
Unfortunately for me, these maps have collars that I want to remove so that I could create a new geotiff file with only the real content of the original geotiff.
I have tried Global Mapper's "Crop collar" command, but that didnt gave me the results I wanted to have.
In my understanding I should use a tool like GDAL warp to complete this task in bulk for all of the images (~106 images in total), but the situation is that I am completely inexperienced using that tool.
I have already tried running gdalwarp -srcnodata 0 -dstalpha input.tif output.tif but that didn't gave me good results as well, therefore I would be glad to receive some help!

Thanks in advance


texts: Creating polygon from exterior points using PostGIS?

texts: How to convert points data to polygon,  so as to form polygon from exterior point as shown above


texts: Adding base layer for WMS layers from WMSCapabilitiesStore

texts:
I'm using the following code to get my WMS layers through GeoExt's WMSCapabilities store and display them in a grid for the users to select and add them to the map panel. but i need to add also a base layer but i couldn't succeed.

 Ext.onReady(function() {
        app = new Ext.Viewport({
            layout: "border",
            items: items
        });
    });
    var map = new OpenLayers.Map(null, {
          projection: new OpenLayers.Projection("EPSG:3857"),
          allOverlays: true,
          maxExtent: new OpenLayers.Bounds(-20037508.34, -20037508.34,
                                 20037508.34, 20037508.34),
          units: 'm',
          numZoomLevels: 10
      });
    items.push({
        xtype: "gx_mappanel",
        ref: "mapPanel",
        region: "center",
        map: map
    });
    controls.push(
        new OpenLayers.Control.Navigation(),
        new OpenLayers.Control.Attribution(),
        new OpenLayers.Control.PanPanel(),
        new OpenLayers.Control.ZoomPanel()
    );
    items.push({
        xtype: "grid",
        ref: "capsGrid", // makes the grid available as app.capsGrid
        title: "Request/Status Map Layers",
        region: "north",
        height: 150,
        viewConfig: {forceFit: true},
        store: new GeoExt.data.WMSCapabilitiesStore({
            url: "data/getcap_1.3.0.xml",
            autoLoad: true
        }),
        columns: [
            {header: "Name", dataIndex: "name", sortable: true},
            {header: "Title", dataIndex: "title", sortable: true},
            {header: "Abstract", dataIndex: "abstract"}
        ],
        bbar: [{
            text: "Add to Map",
            handler: function() {
                app.capsGrid.getSelectionModel().each(function(record) {
                    var clone = record.clone();
                    bbox = clone.get("llbbox");
                    clone.getLayer().mergeNewParams({
                        format: "image/png",
                        transparent: true
                    });
                    var osmLayer = new OpenLayers.Layer.OSM();
                    //app.mapPanel.layers.add(osmLayer);
                    app.mapPanel.layers.add(clone);
                    app.mapPanel.map.zoomToExtent(
                        OpenLayers.Bounds.fromArray(clone.get("llbbox")).transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:3857"))
                    );
                });
            }
        }]
    });

    </script>
</head>
<body>
</body>


texts: DSM quality lower in ArcGIS Pro than PhotoScan

texts:
I am processing aerial photos collected with an automated drone flight into orthomosaics and DSMs, and have noticed the quality of outputs produced in PhotoScan is much higher than those produced in ArcGIS Pro. DSM's produced in ArcGIS Pro have lower surface detail, more noise, and lose significant portions of the target area around the edges of the mosaic dataset.

The process in PhotoScan is fairly automated-- all you have to do is align photos, generate a point cloud, and generate a DSM. ArcGIS Pro provides many more geoprocessing options at every step of the OrthoMapping process, and I've tried nearly all of them. The outputs shown below are the best I've been able to squeeze out of it so far, and were created with the OrthoMapping DEM Wizard using Enhanced Semiglobal Matching, IDW interpolation, no smoothing, and a pixel size equal to 1xGSD (1 cm). Results with other matching methods (Extended Terrain Matching or Semiglobal Matching) and interpolation methods (TIN linear or TIN nearest neighbor) are noticeably poorer. Adding smoothing also doesn't seem to reduce surface noise, discards more real data than is desirable, and results in holes in the finished DSM.

There are many tie points distributed evenly across the entire area (not shown), and the point cloud generated by ArcGIS Pro has 2 million points, which is actually denser than the point cloud produced with PhotoScan. So I don't know why Pro is cutting off parts of the dataset and providing less detail in the areas it keeps.

Are there any parameters that I can change in ArcGIS Pro to improve the quality of the results? If possible, I'd like to stick with Pro-- we already have full ESRI licenses. I'm currently using a 30 day trial of PhotoScan, and would prefer to avoid having to purchase a license for a separate program. On a more practical level, I also appreciate that processing steps in Pro are less of a black box.

Here is the orthomosaic with photo locations (yellow points), flight legs (orange lines), and ground control points (red crosses):

Here is the hillshaded DSM from Photoscan:

Here is the hillshaded DSM from ArcGIS Pro. At this scale, it is easy to see that ArcGIS Pro is discarding data around the edges of the dataset:

Closeup of PhotoScan hillshaded DSM. Just to give you a sense of scale, the lines in the image are individual plow furrows. If you zoom in a little more, detail remains crisp (e.g., tread marks from combine tires are visible):

Closeup of ArcGIS hillshaded DSM. Note the lower level of detail and greater surface noise:


texts:
I'm trying to reclassify a raster based upon two criteria 1. Raster cell on "CBH" is less than or equal to 2 AND "footprint" is not null. Then output raster reclassifies any cell inside this criteria to '3'.

Not sure of input and output raster syntax. What I have so far.

("us_110cbh_1" <= 2)  &  ~IsNull("foot_buff_100").

I used a syntax found here a few weeks ago that included 'inras' and 'outras', but of course I didn't save it and can't find it again.

Thanks.


texts: Maintaining point data attributes when creating polylines?

texts:
I created a model in ModelBuilder to take a set GPS coordinates (containing latitude, longitude, timestamp, user_id) from different users, create a polyline feature class, and then split the polyline to get segments between each point for a given user. Think of it as getting the individual legs of a trip, for an individual user.

Is it possible, when creating the line segments, to maintain the attributes from the vertices of the points in which a segment was created (e.g. coordinates and the timestamp of the beginning and ending vertices)?

The purpose is to eventually calculate the velocity of an individual on a given segment of their trip.

Below is my model in python:

import arcpy
arcpy.CheckOutExtension("Spatial")

#Local variables:
infc = "C:\Path\infc" #input feature class
points2line = "C:\Path\infc_points2line" #polyline feature class
points2line_split = "C:\Path\infc_points2line_split" #line split feature class
points2line_split = "C:\Path\infc_points2line_split"

#Points to line
arcpy.PointsToLine_management(infc, points2line,"user_id", "timestamp", "NO CLOSE")

#Split Line at Vertices
arcpy.SplitLine_management(points2line, points2line_split)

#Add Field for distance conversion
arcpy.AddField_management(points2line_split, "distance_km", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

#Calculate Field, convert Shape_Length from DD to KM
arcpy.Calculate_management(points2line_split, "distance_km", "[Shape_Length] * 111.32", "VB", "")

#Line Density
arcpy.gp.LineDensity_sa(points2line_split, "NONE", infc + "linedensity", "0.3", "0.1", "SQUARE MAP UNITS")
#Cell Size and Search Radius = decimal degree - convert to meters: 1.0(111.32km), 0.1(11.132km), 0.01(1.1132km), 0.001(111.32m), 0.0001(11.132m), 0.00001(1.1132m)


texts: How to use a specific transformation for a coordinate projection with Python

texts:
I try to use Python to project single coordinates from Gauß-Krueger-coordinate system (EPSG:31467) to ETRS89-UTM (EPSG:25832).
I'm using the way Antonio Falciano showed here.

So my script looks like this:

from pyproj import Proj, transform

x_in = 3468511.6749226563
y_in = 5428925.356241324

input_proj = Proj(init='epsg:31467')
output_proj = Proj(init='epsg:25832')

x_out, y_out = transform(input_proj, output_proj, x_in, y_in)
print x_out, y_out

returns 468448.880312 5427193.93989

This works for me, but now I want to add a specific geographic (datum) transfomation (DHDN_To_WGS_1984_4_NTv2). How can I manage this. I would prefere to do this without using a GIS but I cannot find a way so far.

After the tip of mkenedy I found out, that I have to use cs2cs and a .gsb-file.
I think the following line has to put in my code somehow, but I really have no idea how. If I try, it'll cause a SyntaxError: invalid syntax:

cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel + units=m + nadgrids=BETA2007.gsb +to +proj=utm +ellps=GRS80 +zone=32 +nadgrids=@null x_in x_out

Although i have been reading internet pages about cs2cs for hours now, I don't understand how to use cs2cs correctly.
Can someone help me with this?

Now I tried to do an reprojection first like it is done on this site.
It looks like this:

from pyproj import Proj, transform
import pyproj

# GK3 = EPSG: 31467
input_EPSG = '31467'
# UTM32 = EPSG: 25832
output_EPSG = '25832'
gsb_filename = 'BWTA2017.gsb'
x_in = 3468511.6749226563
y_in = 5428925.356241324

gk = pyproj.Proj("+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +nadgrids=" + gsb_filename + "+units=m + no_defs")
x, y = gk(x_in, y_in)
print x, y

returns: 1e+30 1e+30
RuntimeError: latitude or longitude exceeded limits


texts: Will Processing without a Spatial Reference cause Issues?

texts:
I was given a bunch of ASCII files that I needed to convert to raster and then create maps from. For some maps I had to do a simple difference calculation between two rasters (did this in ArcMap).

Well, after I had created dozens of these maps I realized the rasters didn't have a spatial reference. The ASCII files all contain the same spatial extent in lat/lon and same cell size in degrees.  I thought to avoid any potential problems that may be caused by the lack of spatial reference, I would define the projection as WGS 84.

Would the lack of spatial reference have caused any problems when I did the difference calculations? I tried calculating the difference for a couple of sample files with and without a spatial reference and found no difference so I assume I can just add the spatial reference to my difference rasters that don't have a spatial reference instead of redoing them all.

Is this correct?


texts:
I would like to ask you about the following code, it works, but I don't think is very helpful to have elif, elif,elif... always. Then, I have to optimise it. The problem is that I can't access to features in every string of each dir() (or don't know how to do it)

def search(your_search):

dir(iface)
dir(qgis.core)
dir(qgis.gui)
  '''
   i = 0
   list_iface = dir(iface)
   list_core = dir(qgis.core)
   list_gui = dir(qgis.gui)

  #while your_search not in list_iface, ?? oh no...
   '''
    if your_search in dir(iface):
     print '\n'
     print "THE CLASS/FUNCTION/METHOD ", "'",your_search,"'", " IS IN iface API"
    elif your_search in dir(qgis.core):
     print '\n'
     print "THE CLASS/FUNCTION/METHOD ", "'",your_search,"'", " IS IN qgis.core API"
    elif your_search in dir(qgis.gui):
     print '\n'
     print "THE CLASS/FUNCTION/METHOD ", "'",your_search,"'", " IS IN qgis.gui API"
    else:
        print '\n'
        print "IS NOT IN THESE THREE APIs"


texts: Publish GeoTIFF-Layer from Geoserver in Mapproxy

texts:
This is a followup question to this one:

Access WMTS from a GeoTIFF Layer in GeoServer

I'm pursuing approach 1:

Request the GeoServer layer as a WMTS tile set via GeoWebCache (which is built into your GeoServer install) - MapProxy should be able to combine two WMTS layers.

The problem is, I can't get it done. The information found at https://www.geowebcache.org/docs/current/services/wmts.html seems to be largely outdated. I've tried to puzzle the url together myself, alas to no avail.

For reference, the layer in Geoserver is called bel_geo.

Here's my currently nonworking mapproxy configuration:

services:
  demo:
  tms:
    use_grid_names: true
    origin: 'nw'
  kml:
      use_grid_names: true
  wmts:
  wms:
    md:
      title: MapProxy WMS Proxy
      abstract: MapProxy Configuration for BEL

layers:
  - name: combined_layer
    title: geotiff + osm
    sources: [osm_cache, geotiff_cache]
  - name: bel_mapnik
    title: TMS layer from tiles
    sources: [osm_cache]
  - name: bel_geo
    sources: [geotiff_cache]
    title: bel_geo

caches:
  osm_cache:
    grids: [osm_grid]
    sources: [osm_tiles]
    format: image/png
  geotiff_cache:
    grids: [osm_grid]
    sources: [bel_geo_wmts]
    format: image/png

sources:
   osm_tiles:
     grid: GLOBAL_WEBMERCATOR
     type: tile
     transparent: true
     url: http://10.77.70.19/osm_tiles/%(z)s/%(x)s/%(y)s.png
   bel_geo_wmts:
      coverage:
         bbox: [83.25852589091137, 17.64976906844559, 83.38175007681757, 17.728854817689484]
         srs: EPSG:4326
      url: http://10.77.70.19:8080/geoserver/gwc/service/wmts/rest/bel:bel_geo/EPSG:4326/EPSG:4326:%(z)s/%(x)s/%(y)s.png
      type: tile

grids:
    osm_grid:
        srs: EPSG:900913
        origin: nw

globals:

The mapnik layer works just fine, but I'm not getting anywhere with the geotiff.


texts: Èxport DXF from shapefile with layer styles in QGIS

texts:
I have Mastermaps 3 layers with layer styles on QGIS. This is exported with Project -> Import/ Export -> Export project to DXF with the following settings;

However when this is imported to AutoCAD the styles cannot be seen.

When tried with Symbology mode - Symbol layer symbology and Feature symbology the DXF file shows nothing when it is imported to AutoCAD.

Can you please help me to to export shapefile with layer style to DXF?


texts:
From where can I download WKT files for US zip codes?

I need to plot zip code boundary on a map using OpenLayers,


texts: Updating field values of footprint attribute table from an excel sheet

texts:
We have 60 - 70 footprints for which the attributes are to be updated for 25 columns. Some of the columns will have coded values (Eg: Country, Metadata standard, Language). I've slightly modified the code provided in this thread Using UpdateCursor for joined field in ArcPy? (Second one - Code provided by Emil Brundage ). Since i don't know how to make the fields to be read in loop, for 25 fields i've used the same code. As an example i've only listed three fields.

Every column in the footprint is getting updated with the values from an excel sheet w.r.t "Name"column in both footprint & excel sheet. But "Language"column has not been updated and all i am getting is ""Create cursor has failed". "Language""column has coded values (English, Spanish, Portugese etc..). But the other columns with coded values are getting updated without any issues. Any suggestions on how to sort this out ?

import arcpy

footprint = arcpy.GetParameterAsText(0)
footprintFld1 = "Title"
footprintFld2 = "Date"
footprintFld3 = "Language"
footprintIdFld = arcpy.GetParameterAsText(1)

excel = arcpy.GetParameterAsText(2)
excelFld1 = "Title"
excelFld2 = "Date"
excelFld3 = "Language"
excelIdFld = arcpy.GetParameterAsText(3)

valueDi = dict ([(key, val) for key, val in arcpy.da.SearchCursor (excel, [excelIdFld, excelFld1])])
with arcpy.da.UpdateCursor (footprint, [footprintFld1, footprintIdFld]) as cursor:
    for update, key in cursor:
        if not key in valueDi:
            continue
        row = (valueDi [key], key)
        cursor.updateRow(row)

del cursor

valueDi = dict ([(key, val) for key, val in arcpy.da.SearchCursor (excel, [excelIdFld, excelFld2])])
with arcpy.da.UpdateCursor (footprint, [footprintFld2, footprintIdFld]) as cursor:
    for update, key in cursor:
        if not key in valueDi:
            continue
        row = (valueDi [key], key)
        cursor.updateRow(row)

del cursor

valueDi = dict ([(key, val) for key, val in arcpy.da.SearchCursor (excel, [excelIdFld, excelFld3])])
with arcpy.da.UpdateCursor (footprint, [footprintFld3, footprintIdFld]) as cursor:
    for update, key in cursor:
        if not key in valueDi:
            continue
        row = (valueDi [key], key)
        cursor.updateRow(row)

del cursor
...
Runtime error
Traceback (most recent call last):
File "<string>", line 35, in <module>
RuntimeError: create cursor has failed
>>>

35th line in the code is // valueDi = dict ([(key, val) for key, val in arcpy.da.SearchCursor (excel, [excelIdFld, excelFld3])]) //. And, i am running this script as a tool from ArcToolbox.


texts:
I have NIR, green, red edge, red, IR bands at disposal of a an area with a river, I want to discriminate water surfaces, do you have any good ideas on how to do that ? I did a NDWI which is not very convincing. NDVI neither.
Multispectral sensor is Multispec4 from Airinov.

thermal:7.5-13.5 micrometer, sensor is a thermoMAP camera of Sensefly, produced by airinov I think

The NIR/(red+green+NIR) ratio is visible below, river streams are kind of visible but there is lot of noise:


texts: Using Vector Ruggedness Measure (VRM) tool from ArcScripts at latest version of ArcGIS for Desktop?

texts:
I need to use the Vector ruggedness measure (VRM) tool available from http://arcscripts.esri.com/details.asp?dbid=15423. However, I am using ArcGIS 10.2 for Desktop and these tools are made for 9.0, 9.1 and 9.3 only. I have tried to run this tool regardless and get the error:

  Traceback (most recent call last):   File
  "C:\Users\jc221340\Documents\project\Processed
  data\My_tools\Downloaded\ruggedness.py", line 20, in
      gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")
  RuntimeError: Object: Toolbox C:/Program
  Files/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx does not
  exist

  Failed to execute (Ruggedness(VRM)).

How can I use this tool?


texts: How to Unlock feature classes in ArcCatalog that don't show user

texts:
In ArcCatalog there are some feature classes that are locked and don't show a user. After restarting SQL server, the feature classes remain locked. Any ideas as to what might be causing this and how I can unlock them?

Map services should have a user name like gisreader or something attached to them.

Versions: ArcCatalog version 10.1, SQL Server 2008r2 and SDE 10.0.


texts: Size of relation between amount of pixels and file size in GeoTIFF

texts:
I have to calculate upfront from an polygon, how big the GeoTIFF is going to be in bytes. I have a fixed pixel size and know how many pixels there are going to be in the GeoTIFF. The GeoTIFF has only one band.

How is the relation between amount of pixels and file's size in GeoTIFFs?


texts:
I have a geodatabase file with attributes columns that exist in the correct order, however, when I view the attribute table in QGIS it re-sorts that order (possibly due to a previous copy-paste of styles).

Is there any way to reset the attribute table settings for a layer, and let QGIS start afresh displaying the columns in their native order?


texts: I need to change the layer visibility for users in my application. Layers are served by OL 3 by a wms/wfs request. The problem is that when I change the accessibility for users, banning someone to see layers, they still have access to them  as long as they don't log out or close the browser.


texts: QGIS : How to add points from GPS stream data

texts:
I try to place points on QGIS using the GPS plugin but it does not like me.

GPS stream is good, I have the read on the data on port com 66,

QGIS seems to be good with it since it locates correctly where I am, finds latitude, longitude and elevation, satellites and so on:

but saving a point using the add entity on a point vector layer in edit mode returns NULL geometry.

I tried on a shape and spatialite layers.

Should I define a special coordinate system for my layers?

Adding a point manually works fine.


texts: Add point to existing file with mouse click in ArcGIS Pro

texts:
How do I add a point to a point shapefile in ArcGIS Pro 2.0?
All the modify options seemed geared toward polygons. I have looked through the documentation and can't find anything.


texts: Finding maximum distribution between multiple columns

texts:
I have an attribute table that shows how many buildings were built in what decade. The columns show the percentages for every 10 year (1939, 1940, 1950, 1960, 1970, 1980, 1990,2000, 2010) step for every polygon.

This is what it looks like:

Now I need to find out in which of the columns the maximum distribution lies. In the end, I'm supposed to have a choropleth map that shows the main decade it was built in. Does that make any sense?

I've already tried working with the styles but for the categorized settings you can always just choose one column but I need 9. Any tips?


texts: GeoServer returns 401 unauthorized when requesting data via openlayers

texts:
When using the line crossOrigin: 'anonymous',, it returns 401 unauthorized response.

When commenting out that line it works, but I have to be log in to the GeoServer.

Update:
I have changed the Request URL to the below format. Now i get a Deprecation message as shown in the yellow border at the bottom, only in chrome. In Firefox it works.

var geoURL = "//username:password@81.4.../geoserver/web/wms";

wmsSource = new ol.source.ImageWMS({
  url: geoURL,
  params: {
    'LAYERS': "web:streets" , 'TILED': true
  },
  serverType: 'geoserver',
  //crossOrigin: 'anonymous',
});

When accessing that url i get:

  [Deprecation] Subresource requests whose URLs contain embedded credentials (e.g. https://user:pass@host/) are blocked. See https://www.chromestatus.com/feature/5669008342777856 for more details.


texts: QGIS : On layers using my custom function and also having a Rule (Filter) the Vertex Tool is not working

texts:
I am using QGIS 3.16.
To show my problem in details, I have made a very simple test example :
Create New Shapefile Layer
File name :  Test
Geometry type :  Line   OK
Scale :  1:1000
Draw a line     id=1
Layer properties :
.   Rule based
.   Rule (Filter) :  id=1
.   Simple Line :
.   .   Stroke width :
.   .   .   Expression String Builder :
.   .   .   .   Function editor :
.   .   .   .   .   +      New file name :  sf
.   .   .   .   .   Insert my function code (found at the buttom)
.   .   .   .   .   Save and load functions
.   .   .   .   Expression :  sf( 3 )     OK   OK
.   Style :  Save Layer Style :  Test.qml     Save   OK   OK
Notice that it is possible to edit the line with the Vertex Tool
Save Layer Edits
Save the Project File :  Test.qgs
Close QGIS
Restart   Test.qgs
Notice that it is NOT possible to edit the line with the Vertex Tool
Layers properties :
.   Remove the Rule (Filter) :  id=1     OK   OK
Notice that it is now again possible to edit the line with the Vertex Tool
Layer properties :
.   Reinsert the the Rule (Filter) :  id=1
.   Simple Line :
.   .   Stroke width :
.   .   .   Expression String Builder :  Replace  sf( 3 )  with
.   .   .   .   3 / if(@map_scale<1000, @map_scale/1000, 1+ln(@map_scale/1000))
.   .   .   .   which is doing the same calculation as  sf( 3 )     OK   OK
Notice that it is still possible to edit the line with the Vertex Tool
Conclusion :  See title
Anyone having suggestions to solve the problem ?
My function :
import math
from qgis.utils import iface
from qgis.core import *
from qgis.gui import *
@qgsfunction(args='auto', group='Custom')
def sf( x, feature, parent):
.   s = iface.mapCanvas().scale()/1000
.   if s<1 :
.   .   y = x / s
.   else :
.   .   y = x / (1+math.log(s))
.   return y


texts: I have strange results when run Qgis2Threejs: used Google satelite, SRTM file in format .hgt, projections: 3857 and 4326 but didn`t receive 3D model.


texts: Smoothing buildings created from LiDAR to Raster/TIN using ArcGIS Pro?

texts:
I have just started working with LiDAR data. My end result would be a DSM.

To start with creating this I used the tools LAS to Raster and LAS to TIN. With both the results I'm getting an outcome looking like this:

I've made sure it takes the maximum value while converting. Anyone knows how to improve the quality of the created buildings? I am using ArcGIS Pro.


texts: how to get mapnik to render tiles from xml mapfile

texts:
I have Mapnik set up on Ubuntu.  Tirex directs rendering requests from an OpenLayers client browser to Mapnik, which is set up to use the xml file below.

<?xml version="1.0" encoding="utf-8"?>
<Map background-color="#eeebe2" srs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs">

    <Style name="ocean">
        <Rule>
            <PolygonSymbolizer fill="blue" />
        </Rule>
    </Style>
    <Layer name="ocean" status="on" srs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs">
        <StyleName>ocean</StyleName>
        <Datasource>
            <Parameter name="type">shape</Parameter>
            <Parameter name="file">ocean</Parameter>
        </Datasource>
    </Layer>

</Map>

Tiles are generated, but they are all single-color #eeebe2 (from the Map's background-color attribute).  The polygons in the ocean shapefile do not appear.

How can I make the tiles that are generated reflect the polygons in the shapefile?

As an aside, I can run the following simple python script manually, and it outputs the image using the style that is defined in the xml file.  So I am missing some sort of step related to generating the tiles.  Is it a projection thing?  A server setting?

mapxmlfile = 'teststyle.xml'
outputfile = 'outimage.png'
m = mapnik.Map(800, 600)
mapnik.load_map(m, mapxmlfile)
m.zoom_all()
mapnik.render_to_file(m, outputfile)


texts: Drawing lines on Maps (eg. Neighbourhoods or Streets) and then importing into SQL Server?

texts:
I wish to draw some custom boundaries (eg. neighbourhoods or streets or city boundaries).

What software (free or paid) can I use to do this? I read that google maps recently has support for this?
Import this boundary into Sql Server 2008.

Are there any links, tutorials or info anyone can provide, please?


texts: Finding Multistring inside multipolygon using PostGIS?

texts:
If we have below two tables, how can we find data for all multistring(vicRoads_Data table) lying inside multipolygon(patient_data table)?

I tried this

lung_cancer=# select * from patients_geom,vicRoads_Data where ST_Within(patients_geom.geom,vicRoads_Data.geom)

but no luck.

lung_cancer=# \d patients_geom;
                  Table "public.patients_geom"
           Column       |          Type          | Modifiers
    --------------------+------------------------+-----------
     lung_deaths        | real                   |
     heart_deaths       | real                   |
     respiratory_deaths | real                   |
     copd_deaths        | real                   |
     area_code          | integer                |
     area_name          | character(250)         |
     gid                | integer                |
     state_code         | character varying(1)   |
     sla_code06         | character varying(9)   |
     sla_name06         | character varying(50)  |
     sla_5digit         | character varying(5)   |
     geom               | geometry(MultiPolygon) |

    lung_cancer=# \d "vicRoads_Data";
                                          Table "public.vicRoads_Data"
       Column   |           Type            |                           Modifiers
    ------------+---------------------------+---------------------------------------------------------------
     gid        | integer                   | not null default nextval('"vicRoads_Data_gid_seq"'::regclass)
     allveh_amp | character varying(254)    |
     label      | character varying(254)    |
     allveh_pmp | character varying(254)    |
     road       | character varying(254)    |
     trucks     | integer                   |
     allvehs    | integer                   |
     __gid      | integer                   |
     directio_1 | character varying(254)    |
     hf_1       | integer                   |
     year       | integer                   |
     per_trucks | character varying(254)    |
     geom       | geometry(MultiLineString) |
    Indexes:
        "vicRoads_Data_pkey" PRIMARY KEY, btree (gid)
        "vicRoads_Data_geom_idx" gist (geom)


texts: Can't import shapefiles into SQL Server, because of Encoding issues

texts:
I'm trying to import some shapefiles into SQL Server.

At first I had difficulty finding third party tools as I realized that Microsoft has not supported this out of the box. So I found Shape2Sql and Ogr2Ogr.
Then I realized that Shape2Sql is not under active development so that I had to replace some DLL files to make it work.
Then I imported some shapefiles from geonames.com.
Then I found better shapefiles in http://www.openstreetmap.org, downloaded from here
Then I tried to import them using Shaep2Sql, but encoding was not correct.

Then I tried with some different collations to see if the problem goes away or not, and no success.
Then I tried to import it using Ogr2Ogr (version 2.18), but it failed with memory failure error. I searched and found out that this is a reported bug and it's not solved yet.
Then I tried to download older version (before 2.0) and it didn't have ogr2ogr exe file.

And now I'm stuck, spending more than 2 days trying to make it work and the more I search the less I find.

Do we have other tools for importing shape files into SQL Server? Is it possible to do it via C# code through some library? I have no idea cause I'm new to GIS and I don't need much from GIS. Only bounding boxes for a bunch of cities. That's all I want.

Update: I also tried importing dbf files into Excel by simply opening them with Excel, and again Unicode characters are not shown properly:


texts:
I have a vector contour data(polyline). I want to add values to each contour lines. Means I want to digitize my contour lines with proper values. Those values are with me in a paper. I'm adding those (z values) in the "id" column in the attribute table after selecting each line separately. it's taking a lot of time since there are around 20000 lines. Actually I created these contours after digitizing my raster map.(raster to vector conversion after scanning those maps)

Is there any way to add values automatically in between lines, if I select two extreme lines (contour is varying by an interval of 1 unit) or is it possible to create a GUI like that, or can I program like that in ArcMap 10.4?

I'm new to this.

I checked whether is there any pattern of change in "FID" of lines. But I couldn't find a common pattern. I'm adding my contour values to the "id" table.


texts: QGIS 3 Scale Bar Should Print at 1" But Does Not

texts:
In QGIS 3.4, I am using the Print Layout to print some UTM lines with scale bars.  My goal is to make a PDF of a map that is at 1:63360 or 1 inch to the mile.

When I print the PDF, the mile Scale Bar does not measure out to be one full inch.  I have ensured that I am not introducing scaling during the print rendering; that is, I print at 100% scale from AirPrint in the macOS Preview.

Other useful Info

Edit — Root cause found

The problem arose because I copy/pasted the legend from Map A to Map B.  My guess is when pasting the new legend object into Map B preserved the scale from Map A.

I have added three Scale Bars to the layout that appear to be working as expected.

Km - Single Box, colored in yellow
Miles - Line Ticks, in black
Numeric, which shows the scale ratio — 1:63360

Importing the PDF into GIMP and setting the resolution to 300 yields these measurement

Km is 220 px at 300 dpi = 0.73" when printed
Miles is 360 px at 300 dpi = 1.2" when printed  already wrong
But, the ratio of 360 px for 1 mile to 220 px for 1 Km or 360/220 is approximately 1.6. Therefore, the Scale Bar are correct when generated from QGIS.  The measurement errors in pixels is due to me drawing a line in GIMP, so I'm not too worried about that difference.

I've tried the following projections which both yield similar errors of one mile not equal to 1"

epsg:3857
epsg:32611  (WGS 84 / UTM zone 11N).  The UTM data and the rest of the map is from this UTM zone.


texts:
Is there a way to check if data was versioned using "move edits to base" option?

Preferably through ArcGIS itself or using Python, but any way will work.

EDIT

Is there way to see if this box was checked when versioning the data:


texts: Delete points on the map to remain only one for each plane in OpenLayers

texts:
I'm trying to put on the map the route of an ads-b plane, I've managed this so far ...
Poblema is as follows puts all the points since I started to intercept him until the end what is not correct should remain only the last point that shows his exact real time location.
function addMarker(coordinates, data_no) {

            //console.log(coordinates);
            //var marker = new ol.Feature(new ol.geom.Point(coordinates).transform('EPSG:4326', 'EPSG:31700'));
            var marker = new ol.Feature(new ol.geom.Point(coordinates).transform('EPSG:4326', map.getView().getProjection()));
            //console.log(marker);
            var zIndex = 1;

            var style_aircraft = new ol.style.Style({
                                      image: new ol.style.Icon(({
                                        anchor: [0.5, 36],
                                        anchorXUnits: "fraction",
                                        anchorYUnits: "pixels",
                                        opacity: 1,
                                        src: "img/airport.png",
                                        zIndex: zIndex
                                      })),
                                      zIndex: zIndex,
                                      text: new ol.style.Text({
                                        font: 'bold 16px "Open Sans", "Arial Unicode MS", "sans-serif"',
                                        placement: 'center',
                                        fill: new ol.style.Fill({
                                          color: 'black'
                                        })
                                      })
                                    })

            var data = data_no[0]+"/"+data_no[1]+"/"+data_no[2]+"/"+data_no[3];

            style_aircraft.getText().setText(data);
            marker.setStyle(style_aircraft);
            vectorSource.aircraft.addFeature(marker);

          }


texts: Getting extent of composer map in PyQGIS?

texts:
In answering a labeling question: Changing to alternate label if first label does not fit in QGIS? I found myself wondering if there was a way to programmatically get the extent of a composer map so it could be displayed in labels and used in expressions.

So the question arose: Is there a way to find the extent of a specific map on a specific composer programmatically in a way that could be used in the expression engine? I've come up with a solution (I'm posting as an answer below) But I'd be happy if there was a more efficient solution. If you have a better one let me know.


texts: Extract morphologic characteristics of mountains from DTM in QGIS

texts:
I'm looking for some algorithm/script that allows the automatic extraction of peaks from DTMs into a point shapefile (including their morphologic characteristics such as key col, prominence, and isolation, as indicated in the image). I've tried r.geomorphon but it returns only a raster with morphologic characteristics of the entire DTM.
I've always done it manually, but it's really time-consuming for thousands of features.


texts:
Is there a way of showing the progress when adding a layer though QgsVectorLayer(url, name, "ogr") in pyqgis? Some of the layers I am adding are quite large, and take some time to load.

Excample:

    url = "http://wfs.geonorge.no/skwms1/wfs.tilgjengelighettettsted?service=WFS&request=GetFeature&version=2.0.0&srsName=urn:ogc:def:crs:EPSG::4258&typeNames=app:TettstedHCparkering&"

    new_layer = QgsVectorLayer(url, "LayerName", "ogr") #this part might take long time

    QgsMapLayerRegistry.instance().addMapLayer(new_layer)


texts: HEC GeoHMS generates wrong catchment area

texts:
I am trying to define a catchment area automatically with HEC Geo HMS (ArcGIS version 10.1).
My input-data is a DEM-file (raster 1m). After doing all steps in HEC Geo HMS, the catchment area (see black line in the picture) is generated.

Comparing this catchment area with the contours, I think there is a mistake. In my opinion the northern boundary of the catchment area should reach the highest elevation point of my DEM (which is marked by an orange contour). Can anybody imagine, where the error in the procedure is?


texts: Batch conversion of XY coordinates to latitude and longitude

texts:
I'm working on a school project. I have this dataset from Sacramento, California

I want to convert the XY coordinates to latitude and longitude, in order to import it into tableau and use its mapping features.

I also have the address, if that makes it easier.

I have about 10,000 records, so I need a powerful tool to perform this conversion.


texts: Why is mosaic of rasters big, black rectangle in ArcGIS Desktop?

texts: I downloaded NED rasters from USGS, and I now want to mosaic all of the rasters into one. I used the Mosaic to New Raster tool, but the output is a big, black rectangle. Am I missing something? I feel like I might be making a rookie mistake.


texts: This is a silly question, but one I've been trying to figure out for a while.  I have two Aster DEM tiles, one next to the other.  How do I ensure uniform colour for both?


texts:
I am trying to convert SQL Shapes to ESRI Shapefiles using the instructions in this blog post from a few years back and OSGeo4W.

https://alastaira.wordpress.com/2011/09/29/exporting-spatial-data-from-sql-server-to-esri-shapefile/

I'm most of the way there, but I'm losing a shape here and there on account of some of them being GeometryCollection shape types.  Is it possible to convert GeometryCollections to Polygons on the fly and get all of these shapes out in the same shapefile at once?  Or should I just extract them separately and match them up later?

Not sure it's necessary, but this is the command I'm using in the OSGeo4W Shell

ogr2ogr -f "ESRI Shapefile" "C:\FolderName\AllShapes.shp"
"MSSQL:server=MyServer;database=ShapeDatabase;table=.dbo.Shapes;
trusted_connection=yes;" -sql "Select ShapeID, ShapeTypeID as ShapeType,
ShapeData From .dbo.Shapes Where ShapeData.STGeometryType() In
('Polygon','Multipolygon','geometrycollection')"

That command works if I dump out 'geometrycollection', but, as I mentioned, I'm then missing a bunch of shapes.


texts:
I want to find geolocation with GSM Cell information. I need tower information to find locations and calculate triangulation. I tried OpenCellID service but I think it's not up to date. I can't find the cell towers near me on the OpenCellID map.

Can anybody suggest an up to date free or paid service? I can also accept any other solutions?


texts:
When the tooltip is activated on a geographical layer it pops up with the cursor on the top-left corner of the label thus covering part of it ...

Is there any way to fine-tune this behaviour (without rewriting a select tool in python...)?