texts: Filtering polygon using Cesium for WMS 1.3.0

texts:
I am trying to display only the portions of my data that fall within the polygon. I am using WebMapServiceImageryProvider to get the tiled imagery.
As Cesium just supports WGS84, I tried to use both EPSG:4326 and CRS:84 with different axis ordering. But it was not successful.
This is the way the polygon is represented:
<ogc:Intersects>
  <ogc:PropertyName>gsmlp:shape</ogc:PropertyName>
    <gml:MultiPolygon srsName="urn:ogc:def:crs=EPSG:4326">
      <gml:polygonMember>
        <gml:Polygon srsName="EPSG:4326">
          <gml:outerBoundaryIs>
            <gml:LinearRing>
            <gml:coordinates decimal="." cs="," ts=" ">-16.647690384169568,137.98786059183666 -26.030291233429402,140.9468908103431 -31.7376727729486,128.9819274324566 -32.08022707854407,115.90510193754302 -14.63565118294481,128.98253730075106
            </gml:coordinates>
          </gml:LinearRing>
        </gml:outerBoundaryIs>
      </gml:Polygon>
    </gml:polygonMember>
  </gml:MultiPolygon>
</ogc:Intersects>

My code in “Cesium Sandcastle” is as follow:
var viewer = new Cesium.Viewer("cesiumContainer");
var rectangle =Cesium.Rectangle.fromRadians(2.021479979302771, -0.5587749816864094, 2.460402722308602,  -0.25865874344230144);
var imageryLayers = viewer.imageryLayers;

imageryLayers.addImageryProvider(
  new Cesium.WebMapServiceImageryProvider({
    url:"https://nationalmap.gov.au/proxy/http://geoserver.nationalmap.nicta.com.au/geotopo_250k/ows",
    layers: "Hydrography:bores",
    rectangle: rectangle,
    parameters: {
      transparent: true,
      format: "image/png",
      VERSION: "1.3.0",
      service: 'WMS',
      tiled : true,
      request: 'GetMap',
      srs: 'EPSG:4326',
      sld_body: _sld_body
    },
  }));

viewer.camera.setView({
  destination: Cesium.Rectangle.fromDegrees(
    114.591,
    -45.837,
    148.97,
    -5.73
  ),
});

I used the following code to calculate the polygon’s coordinate and its surrounding rectangle:
const coords = cartesian3
     .map(cart => Ellipsoid.WGS84.cartesianToCartographic(<Cartesian3>cart))
      .map(latLon => [latLon.latitude * 180/Math.PI , latLon.longitude * 180/Math.PI]);

const polygonHierarchy = new PolygonHierarchy(<Cartesian3[]>cartesian3);
const rectangle=PolygonGeometry.computeRectangle({polygonHierarchy:polygonHierarchy , ellipsoid:Ellipsoid.WGS84  });


texts:
I recently used osm2po (v 5.0.0) on SF with these parameters:
java -Xmx8g -jar osm2po-core.jar cmd=tjspg prefix=sf tileSize=x california-latest.osm.pbf

and everything looks great in the whole city except for some streets around Mission and 2nd. I've added a picture for illustration, streets that are not colored are not in the sf_2po_4pgr table. Is there something I'm missing in terms of params - maybe different tileSize?

The only change I've done in the config is added a vertex writer:
postp.1.class = de.cm.osm2po.plugins.PgVertexWriter


texts: Different watersheds for same area because of "burned streams"

texts:
Im trying to analyse where water is flowing which rains down to my project area. This is the DEM I have (my project area is bordered blue, the red marks are my pour points):

I learned from another thread (Could bridges prevent ArcMap from creating correct watersheds?) that culverts could disturb the creating of watersheds. So I used the "DEM Reconditioning" tool of ArcHydro to burn the streams:

Then I used this DEM (called it "agreedem") and the following workflow to generate my watersheds (snap distance set to 3 m):

This was my first result:

Since it is my main aim to analyze what I have to do to keep the water inside the center of my area I started building small dams using the "Build Walls" tool of ArcHydro. After that I rerun my watershed workflow (several times with different dams, this was the best result I could get):

I was pretty happy with the result. Then I asked myself: "Why should I use the "burn the streams" thing in this case?". There are no streams with culverts in this area (though I have the situation with the culverts in other areas I'm looking at) so there should be no reason to do so.
I deceided to rerun the watershed workflow for my original DEM (called "dtm_orig"):

As you can see, the watersheds are different from the ones I generated with the "agreedem" DEM. So I asked myself how there could be such a great difference in the watersheds between the ones created basing on the "agreedem" and the ones based on "dtm_orig". I couldn't imagine that bigger and deeper ditches / streams (whats the only difference between "agreedem" and "dtm_orig") lead to this results and analyzed all the steps of my "create watershed" workflow.

I then realized that there are huge differences in the DEMs which are generated in the first step by using the "Fill" tool between the "agreedem"...

...and the "dtm_orig" DEM:

This are the values of one of the points where one could see the difference:

Why is there such a big difference between the both filled DEMs?

As suggested by Martin I ran the Sink Evaluation (on the dtm_orig DEM) today. To be honest I don't know how to interpret the results.

This is the "SinkPoly" Layer:

Same with added "SinkDA" Layer:


texts: Reading NARR reanalysis data with raster package makes R crash

texts:
I'm trying to analyse data from NARR (NCEP North American Regional Reanalysis) and every time I try my R session crash.

Here is what I do as a reproducible example:

library(raster)

the_NetCDF <- tempfile(fileext = ".nc")

download.file("ftp://ftp.cdc.noaa.gov/Datasets/NARR/Dailies/monolevel/csnow.2019.nc", the_NetCDF)

the_ras <- brick(the_NetCDF)

And here is what I get every time I run the code:

I also tried to open the file in QGIS with also an error.

here is my session info if it can help:

R version 3.6.2 (2019-12-12)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)

Matrix products: default

locale:
[1] LC_COLLATE=English_Canada.1252  LC_CTYPE=English_Canada.1252
[3] LC_MONETARY=English_Canada.1252 LC_NUMERIC=C
[5] LC_TIME=English_Canada.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] ncdf4_1.17    raster_3.0-12 sp_1.3-2

loaded via a namespace (and not attached):
[1] compiler_3.6.2   tools_3.6.2      Rcpp_1.0.3       codetools_0.2-16 grid_3.6.2
[6] lattice_0.20-38


texts: Does anyone know the SRID for an equal area projection that is suitable for Finland?


texts: How to save a joined spatialite table when each table contains PKUID fields?

texts:
QGIS 1.8.0 --
QSpatialite 6.0.7 --
Windows XP

I have spatialite tables joined together that I need to save as a single table, but an error-message says that QGIS can't save the table because there are two PKUID field columns.

Spatialite_GUI can delete most spatialite table columns, when I need to do that to manually adjust sqlite column names, but my Qspatialite or QGIS always generates PKUID columns for every table it makes, and that column is always locked (I suppose as a primary key column)  The PKUID column can't be deleted by Spatialite_GUI since it's locked, and so I get stuck and have to export everything to CSV files. That's the only workaround I know of, to save the table out as CSV file, manually delete the column from CSV and reload it, then join and "save as".

Is there an easier way to save joined tables in QGIS (that have duplicate columns) or locked columns, rather than manually deleting or renaming identically named columns beforehand with CSV files.  This is a chore since sqlite db aren't designed to have their field columns renamed or deleted, and QGIS doesn't automatically rename duplicated column names.

What's the easiest way to save a joined table with duplicate column names in QGIS?  Thanks in advance!


texts: Default password for geoserver 2.1.4 not working

texts:
I have installed geoserver 2.1.4 on Ubuntu 10.04 i have jetty running on port 90.

I get to the login screen. I introduce user: admin password: geoserver and i am getting that user or passowrd are not correct!!

It is supposed that those are default user and pass. Can i change it on a file or something?

Maybe geoserver is not properly installed?

Regards,


texts: Failure to Access DBMS Server Error - ArcGIS/Postgres

texts:
If connected to an SDE enabled Postgres database via ArcGIS (10.1) for a long period of time, maybe an hour +, I will often receive an error from ArcGIS stating "Failure to access the DBMS Server[my schema name.sde.GDB_items]. I can't pin down any specific event causing the error, and my fix is to simply restart ArcGIS. Disconnecting and re-connecting the database connection without an ArcGIS restart results in the same error message.

Its a minor inconvenience, but one that would save myself and others daily grief if corrected. Any ideas?


texts: Metadata for Bing Maps Aerial in ArcGIS Desktop?

texts:
Is there a way to dynamically pull up metadata on the Bing Maps Aerial web mapping service in ArcGIS?  I am working at a 1:245 scale across various locations in OR and CA and need to know if Bing's coordinate system is accurate to within 0.5 meters.  Additionally, I need to know the resolution.  ESRI's description is pretty vague:

  This web map contains the Bing Maps aerial imagery web mapping
  service, which offers worldwide orthographic aerial and satellite
  imagery. Coverage varies by region, with the most detailed coverage in
  the USA and United Kingdom. Coverage in different areas within a
  country also varies in detail based on the availability of imagery for
  that region. Bing Maps is continuously adding imagery in new areas and
  updating coverage in areas of existing coverage.


texts:
I have an .las file, intended to be a DTM (bare earth) that I generated from a photogrammetry software package that has come out just a bit too lumpy. I now have the license for that software tied up for a long time on the next one, and am looking for tools to smooth out this DTM in ArcGIS. I have harnessed my education as best as I could, from generating a meaned grid with LAS to Raster to using Focal Mean and Focal Minimum, generating a raster to contours and back to raster.

It seems that no matter what I do, how large I make the search area or grid spacing, I can't get rid of the lumps. The area I'm working on is very mountainous, and I believe that to be part of the problem.

Does anyone have any really good tips or tools that might help me out here? I have access to an ArcGIS 10.2 advanced license, Spatial Analyst and 3D Analyst. I also have access to Global Mapper.


texts: Obtaining Elevation for single point from ArcGIS REST API?

texts:
I have spent quite a while trying to get the ESRI API working. Basically we would like to pass a lat, long to a service and obtain the elevation back.

Our previous service stopped working - we were using:
https://utility.arcgis.com/usrsvcs/servers/55efcc79e42941e8a98831f987ca67c0/rest/services/WorldElevation/TopoBathy/ImageServer/identify?f=json&geometry={x:20,y:20,spatialReference:{wkid:4326}}&geometryType=esriGeometryPoint&f=json&returnCatalogItems=false

this service now appears to have gone.

I have looked at the ESRI documentation but alas it is really poor
https://developers.arcgis.com/rest/elevation/api-reference/summarize-elevation.htm
should work but just submits a job
Using https://elevation.arcgis.com/arcgis/rest/services/WorldElevation/TopoBathy/ImageServer
Displays in ArcGIS Pro/ArcMap but you can only identify, looking at fiddler the response/request isn't easily reproducible.

Alternative options are Open Elevation which is the opposite of ESRI documentation it couldnt be simpler but alas the baseline dataset isnt high enough quality:
https://api.open-elevation.com/api/v1/lookup?locations=41.161758,-8.583933

Does anyone have any working GET/POST requests to get Elevation details for point data from ESRI?


texts: Error 00368 using arcpy.SelectLayerByLocation_management

texts:
I'm new to programming and I'm having trouble with my script:

#Set join parameter variables
inFeatures = arcpy.GetParameterAsText(0)
inField = arcpy.GetParameterAsText(1)
joinTable = arcpy.GetParameterAsText(2)
joinField = arcpy.GetParameterAsText(3)
#Join road feature class to buffer distance table
arcpy.JoinField_management(inFeatures, inField, joinTable, joinField)
#SCRIPT WORKS TO THIS POINT
#Select roads that intersect the stdsm feature class
stdsM = arcpy.GetParameterAsText(4)
arcpy.SelectLayerByLocation_management(inFeatures, "INTERSECT", stdsM)
arcpy.RemoveJoin_management(inFeatures)

Running it as a script tool in ArcGIS I get this message:

Could someone help me with this?


texts:
I am new in the arcgis runtime. I am trying to develop an android application.

Steps that i follwed

Created an account in arcgis.com .
Then i created one Feature Layer (hosted) from the shape files zip .
Then i share this to everyone.
I copied the url given at the below. url is . Please see the output for the url.
Then i try to implemet means to show this in the map on my android application And trying to do some feature-layer-query.Refference is in the GITHUB
I am trying to replace the service_url to my url. I am doing the search operation and getting the result from the online then trying to zoom that extent. But is noe even showing my map too.
Now i am getting problem. it showing the error in

  05-10 16:31:23.516 20030-20030/com.esri.arcgisruntime.samples.featurelayerquery E/Feature layer query: Feature search failed for: kulathoor. Error=com.esri.arcgisruntime.ArcGISRuntimeException: Invalid JSON

Please see the below image for ref

My code try

@Override
    protected void onCreate(Bundle savedInstanceState) {
        super.onCreate(savedInstanceState);
        setContentView(R.layout.activity_main);

        // inflate MapView from layout
        mMapView = (MapView) findViewById(R.id.mapView);

        // create a map with the topographic basemap
        final ArcGISMap map = new ArcGISMap(Basemap.createTopographic());
        // set the map to be displayed in the mapview
        mMapView.setMap(map);

        // create feature layer with its service feature table
        // create the service feature table
        //mServiceFeatureTable = new ServiceFeatureTable(getResources().getString(R.string.sample_service_url));
        mServiceFeatureTable = new ServiceFeatureTable("https://services7.arcgis.com/7FyZZrSIYfiWYztL/arcgis/rest/services/MyGisFileLayerTest/MapServer");
        // create the feature layer using the service feature table
        mFeaturelayer = new FeatureLayer(mServiceFeatureTable);
        mFeaturelayer.setOpacity(1.0f);
        //override the renderer
        SimpleLineSymbol lineSymbol= new SimpleLineSymbol(SimpleLineSymbol.Style.SOLID, Color.RED, 1);
        SimpleFillSymbol fillSymbol = new SimpleFillSymbol(SimpleFillSymbol.Style.SOLID, Color.YELLOW, lineSymbol);
        mFeaturelayer.setRenderer(new SimpleRenderer(fillSymbol));

        // add the layer to the map
        map.getOperationalLayers().add(mFeaturelayer);

        // zoom to a view point of the USA
        mMapView.setViewpointCenterAsync(new Point(-11000000, 5000000, SpatialReferences.getWebMercator()), 100000000);
    }

Can you find out where i am doing mistake . Why the may is not shwing in my app. Below is the screen shot of my application.


texts: Undefined reference to mapnik::hit_grid from C++ compile?

texts:
I am trying to compile a c++ mapnik program with code in question being:

mapnik::grid my_grid(256, 256, "id");

and the command to compile being:

g++ grid_gen.cpp -std=c++11 -lmapnik -lcairo -lfreetype -lpng -licuuc -lz -ldl -Wall -o grid_gen

which ends up throwing this error:

  grid_gen.cpp:(.text+0x744): undefined reference to
  `mapnik::hit_grid<mapnik::gray32s_t>::hit_grid(unsigned long, unsigned
  long, std::__cxx11::basic_string<char, std::char_traits<char>,
  std::allocator<char> > const&)' collect2: error: ld returned 1 exit
  status

I am on Ubuntu 16.04.2 with g++ 5.4.0.

How do I compile this without throwing any error?

If you are curious, I am trying to generate utf8 grid for the map tiles from c++.


texts: How to represent Typhoon(s) wind speed using Quantum GIS / PostGIS?

texts:
I am not completly new to GIS (was taught ArcGIS at Uni) but I would still consider myself (especially compared to all these experts here)  a newbie in QGIS or PostGIS.

I want to display the windspeed of typhoons with the help of buffers along their tracks. I used PostGIS to combine the points of the typhoons to Lines ( like in this post).
Anyhow, if I buffer the lines, buffer segments from the same typhoon overlap or leave gaps. And I need to change that (Here is a link to a picture showing my problem)

My idea was to create a single line for each typhoon, buffer the line and then split the buffer with the help of my point layer. I need the information on windspeed broken down to the points rather than just an average for the complete typhoon track.

Is there a way to split polygons with the help of point layers? And if not, a way to somehow interpolate the buffer segments that I have right now to close gaps and cut overlapping parts? Doing this manually is not an option at I have to work with typhoon data from 1945-2013.

EDIT:
So... if I buffer the lines in Quantum I kind of"solved" the problem of gaps. So I "just" need to interpolate/cut the overlapping parts at the beginning/end of each segment (see picture: http://postimg.org/image/maca8nrc7/ )


texts: Raster layer not visible after changing color rendering in Qgis

texts:
I created a new raster and it appears in the layer list, its visible on the screen and the values seems to be correct.
When I change the "band rendering" from "Singleband gray" to Singleband pseudocolor" the raster seems to disappear from the screen. The values are still there, possible to find by clicking with the "identify features" arrow.


texts: How to transfer saved colors from computer to computer?

texts: I have a ton of colors saved on the color palette and I was wondering how to transfer this color palette from a machine using QGIS 1.10 to a machine using QGIS 1.12.


texts: Postgis query to convert UTM to Lat long

texts:
I want Postgres/PostGIS queries to convert UTM coordinates from UTM 43N projection to Decimal Degrees 32643 projection.

select ST_X(the_geom), ST_X(the_geom) from table


texts: How to get latest picture of my city from sentinel 2?

texts:
I was trying to get latest picture of my city from Sentinel 2, but had no success.
I have registered on Copernicus website, found my city and drew a square around it.
Then I clicked the "search" icon, but all pictures it turned up seem to be "broken" in different ways:

partial (triangle-like, or trapezoid-like)
totaly black
totaly white
white with acidy-green (not grass) pieces in it.

I thought that I had choosen a wrong sensor, but in a downloaded archive I found JPEG2000 files for different bands - including bands 2,3,4 which I need.
My question is, how to download a google maps-like RGB picture from Sentinel 2.
I have heard that I can get 300m/pixel picture resolution. Is it right?


texts: Converting a CSV with a WKT column to a Shapefile

texts:
Using ogr2ogr in Python, I am trying to convert a CSV to a shapefile. One column in the CSV, named "Polygon", contains WKT like this: POLYGON((long lat, long lat, long lat, etc.)) Currently, I can make a polygon shapefile with the correct projection, but the geometry is empty. How can I modify my ogr2ogr arguments to correctly create the geometry using the WKT in each row? Currently I have something like this:

ogr2ogr.main(["", "-a_srs", "EPSG:4326", "-f", "ESRI Shapefile", "output.shp", "input.csv", "-nlt", "POLYGON"])


texts: Creating database for climate data with MySQL

texts:
I'm working on a project using climate data and MySQL. Before I start to work, what's your advice on the best way to complete this task.

I plan on using the following workflow:

  Analyse the climate-data, creating a mysql database-schema that could hold the data. (all tables should use UTF-8 as charset).
  Download all the raw-data from worldclim.org, using the 30 arc-seconds (~1 km) resolution.
  Import the data. Providing a dump

I have already downloaded the data for Min.Temperature, and it looks like this when I open one of the folders

The data is in .adf format: how can I insert that format into MySQL?

Also, what do you think is the best way to organize the database in MySQL for that kind of data?

Thank you.

UPDATE

I explored and found out that is option in AutoCAD Map 3D (FDO) where I can insert shapefiles in MySQL. Maybe then is the best way to convert coverages to shapefiles and then in Autocad Map 3D to MySQL?

I need to use MySQL. So I'll download . bil files. gissolved said I can save it in a blob field in my database and KevinMayall said I should convert it to ASCII and then to MySQL. These are raster files. What is the best way to import .bil files to database then?


texts: Drawing outline of points in ArcGIS for Desktop

texts: Simple question, can Arc take a shapefile of points, calculate which ones are at the "edges" and draw a polygon or polyline around it. I know there are masks and bounding geometry and all of that, but it just needs to draw a line around and can intersect the exterior points.


texts: Print Mapbox XYZ layer using Mapfish/Geoserver

texts:
I'm having trouble getting a mapbox XYZ layer to work with the Geoserver Print plugin - I just get the following error:

spec.layers[0].baseURL has an invalid value

and a bit further down:

Illegal character in path at index 79

My URL is like as so (some parts redacted for privacy):
https://api.mapbox.com/styles/v1/*********/*************************/tiles/256/{z}/{x}/{y}?access_token=ACCESS_TOKEN

I assumed that the illegal character might be the curly brace, but of course this is required in the URL.

If I remove this layer, the print works correctly.

The layer is added to the spec using this object:

 {
    type: "XYZ",
    baseURL: XYZ, //the url shown above
    extension: "png",
    maxExtent: [0, 0, 700000, 1300000],
    tileSize: 256,
    resolutions: [resolutions], //current map resolution, i.e [40]
    path_format: '${z}/${x}/${y}'
  }

I have also added the host to the whitelist in the yaml config:

hosts:
- !dnsMatch
  host: api.mapbox.com

Geoserver Version: 2.13.0

Has anyone had any success getting this to work?


texts: Adding Spatial Index for Geometric Network?

texts:
I am trying to create Spatial Indexes on all the feature classes in one of my databases.  It has worked fine on most feature classes so far, but when I started on feature classes in a feature dataset suddenly I was getting errors about locks.

Runtime error  Traceback (most recent call last):
  File "<string>", line 2, in <module>   File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\management.py", line 6448,
  in AddSpatialIndex
raise e ExecuteError:
  ERROR 000464: Cannot get exclusive schema lock.
  Either being edited or in use by another application. Failed to execute (AddSpatialIndex).

I tested on a couple of other feature datasets, one worked, one gave the same error.  Both where it failed contained geometric networks.  In a test I deleted the geometric network (leaving all the same feature classes) and the Add Spatial Index tool worked as expected.

Why would the Add Spatial Index fail on a geometric network?  I cannot find anything in the Geometric Network and Spatial Index (and related) docs to indicate why this would happen.

Note that I am able to add a spatial index through the right-click feature class properties.

I am using ArcGIS Desktop 10.5 with Advanced license.  These feature classes are stored in a local file geodatabase.

In response to comment by @Alex Tereshenkov I have checked there are no locks.  I also created a new File Geodatabase and copy/pasted a Feature Dataset that I hadn't used today.  This one has feature classes both part of and separate of the geometric network.  I tested using tool from ArcToolbox to rule out arcpy being difficult.

Any feature class within the geometric network fails

Any feature class external to the network succeeds (all are within the same feature dataset)

As recommended by @Alex Tereshenkov, I tried everything completely new - I created a new FGDB, Feature Dataset, 2 Feature Classes (1x point, 1x line), and added half-dozen features.

I then closed out and created new Geometric Network with those 2 feature classes.

Add Spatial Join still fails on both feature classes.  I have tried with my lines both as Complex (as in my original) and Simple edges, but again it fails.


texts: Import CAD Annotation into an existing Geodatabase Annotation Feature Class?

texts:
I'm using ArcGIS Desktop 10.5.1, and I need to automate (with Python using arcpy) importing the annotation of many DWG files into the same file geodatabase annotation feature class. I can't use the ImportCADAnnotation tool because it always wants to create a new annotation class, and I don't want to overwrite the annotation feature class each time (or create new ones each time).

I've been attempting to use SearchCursor on the CAD annotation and InsertCursor on an empty GDB annotation feature class, and I can populate values for most of the fields (I've made the annotation shape into a Polygon which the GDB seems to need), but I don't know how to populate the Element field of the GDB annotation (It stores a BLOB field of the symbology and text geometry). Is this possible?


texts: QGIS calculate distance of point along a line

texts:
I have a polyline (movement path) and points recorded along the line. I would like to calculate for each point the distance along the line in Quantum GIS. I hope for a result which will add a new column to the point attribute table indicating the distance along which it is on the line. I'd appreciate any thoughts.

Mike.


texts: How to create shapefiles with z coordinate with QGIS

texts: I'm working with QGIS, and I want to create a shapefile with the elevation coordinate, I can't find something similar to what I was used to do in ArcGIS, can you help me?


texts: Creating script tool that will create copy of feature class and offset it by given distance using ArcPy?

texts:
I want to duplicate a polygon feature class and offset all of the polygons by about 10 feet in both the x and y directions. I asked if there was any way to do this last week, and I was informed that I would most likely need to make my own python script using arcpy. I made my own script using arcpy, but it isn't working:

import arcpy
from arcpy import env
import os

env.overwriteOutput = True

# Get arguments:
#   Input polygon feature class
#   Output polygon feature class
#
inputFeatureClass = arcpy.GetParameterAsText(0)
outputFeatureClass = arcpy.GetParameterAsText(1)
xShift = arcpy.GetParameterAsText(2)
yShift = arcpy.GetParameterAsText(3)

shapeName = arcpy.Describe(inputFeatureClass).shapeFieldName

# Create the output feature class with the same fields and spatial reference as the input feature class
arcpy.CreateFeatureclass_management(os.path.dirname(outputFeatureClass), os.path.basename(outputFeatureClass), "POLYGON", inputFeatureClass, "", "", inputFeatureClass)

# Create a search cursor to iterate through each row of the input feature class
inrows = arcpy.SearchCursor(inputFeatureClass)
# Create an insert cursor to insert rows into the output feature class
outrows = arcpy.InsertCursor(outputFeatureClass)

# Create empty Point and Array objects
pntArray = arcpy.Array()
partArray = arcpy.Array()

# Loop through each row/feature
for row in inrows:
    # Create the geometry object
    feature = row.getValue(shapeName)

    partnum = 0

    # Count the total number of points in the current multipart feature
    partcount = feature.partCount

    while partnum < partcount:
        part = feature.getPart(partnum)
        pnt = part.next()
        pntcount = 0

        # Enter while loop for each vertex
        #
        while pnt:

            pnt = part.next()
            shiftedPoint = arcpy.Point()
            try:
                shiftedPoint.ID = pnt.ID
                shiftedPoint.X = pnt.X + float(xShift)
                shiftedPoint.Y = pnt.Y + float(yShift)
            except AttributeError:
                continue
            #shiftedPoint = arcpy.Point(float(pnt.X) + float(xShift), float(pnt.Y) + float(yShift))
            pntArray.add(shiftedPoint)
            pntcount += 1

            # If pnt is null, either the part is finished or there is an
            #   interior ring
            #
            if not pnt:
                pnt = part.next()
                if pnt:
                    arcpy.AddMessage("Interior Ring:")
        # Create a polygon using the array of points
        polygon = arcpy.Polygon(pntArray)

        # Empty the array for the next run through the loop
        pntArray.removeAll()

        # Add the polygons (or 'parts') to an array. This is necessary for multipart features, or those with holes cut in them
        partArray.add(polygon)
        arcpy.AddMessage("Added a polygon to the partArray!")
        partnum += 1

    # Create a new row object that will be inserted into the ouput feature class. Set newRow = row so that it has the same attributes
    # Set newRow.shape = partArray so that the only thing different about this new feature is that its geometry is different (shifted)
    newRow = row
    newRow.shape = partArray

    outrows.insertRow(newRow)

    # Empy the array for the next run through the loop
    partArray.removeAll()

del inrows, outrows

I keep getting this error on line 70

<type 'exceptions.ValueError'>: Array: Add input not point nor array object

I can't figure out why it's giving me this error, since I defined the input as an array.

Does anyone know why I'm getting this error?


texts:
is there any way to convert the scripts which is produced in QGIS as python scripts from processing modeler to codes using osgeo library codes?

i want some codes that can be used in QGIS algorithm in OSGEO library


texts: Refreshing values in attribute table automatically when geometry changes?

texts:
I can't find a solution for auto-refreshing attribute table data.

For instance, how can the data related to features of a Shapefile (e.g., the area) be updated if the feature geometry is changed.


texts: Problem with R Package Rgdal and writing MapInfo File

texts:
I have a problem using "writeOGR" with the R Package Rgdal.
The strange thing is, that this worked just a few weeks ago, but when I tried to start the R script today, it didn't work.

The Error is: "Error in writeOGR(dat, dsn = "Zones.tab", layer = "Zones", driver = "MapInfo File") : Failed to create feature"

The MapInfo Files are stored, but they are empty (just the column Names are there).

I wrote this, so you guys can have a look. Same Error when I start the code:

library(rgdal)
setwd("G:/Löschen")

(dat <- data.frame(Number=c(1,2,3,4,5,6),
                    East=c(371008,371124,371346,371276,371051,371344),
                    North=c(5720763,5720567,5720761,5720493,5720506,5720490),
                    BW=c(320,350,300,380,275,390)))

coordinates(dat) <- ~East+North

proj4string(dat) <- "+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

writeOGR(dat, dsn = "Zones.tab", layer = "Zones", driver = "MapInfo File")


texts:
I have a line layer which is cut at the level of the point layer,

I want it to be connected at the nearest point while saving the forme ( just adding the missing part the that line without moving it) as showing in the next picture

My idea is to calculate the distance between each extremity of the line with the points and take the nearest point then draw a line between them..But the drowing come out bizzarly.

I was wondering ,Is there a simple way of doing this or maybe could anyone help me with the following code?

from qgis.core import *
from qgis.gui import *
import processing
layer = iface.activeLayer()
feat1=layer.selectedFeatures()

for lyrr4 in QgsMapLayerRegistry.instance().mapLayers().values():
  if lyrr4.name().startswith( 'Boitier_B1'):
     ch = lyrr4
     pr= ch.dataProvider()
chs= ch.getFeatures()
for f in feat1:
 min1=10000.55
 min2=10000.55
 attrs = f.attributes()

 idc = f["ID"]
 x1=f["x1"]
 y1=f["y1"]
 x2=f["x2"]
 y2=f["y2"]
 print "id_cab"
 print idc
 print "x1"
 print x1
 print "y1"
 print x1
 print"x2"
 print x2
 print "y2"
 print y2
 if  "%.2f" % x1 != "%.4f" %xc and "%.2f" %x2 != "%.4f" %xc:
     point1 = QgsPoint(x1,y1)
     point2 = QgsPoint(x2,y2)
     for ch in chs:
         id = ch["ID"]
         xc= ch["xcoord"]
         yc= ch["ycoord"]
         print "id_ch"
         print id
         print "xc"
         print xc
         print "yc"
         print yc
         pointc = QgsPoint(xc,yc)
         geom = f.geometry()
         distance = QgsDistanceArea()
         m1 = distance.measureLine(pointc, point1)
         m2 = distance.measureLine(pointc, point2)
         print "m1"
         print m1
         print "m2"
         print m2
         fet = QgsFeature()
         if m1<min1:
             min1=m1
             pts=pointc
         if m2<min2:
             min2=m2
             idmin2=id
             pts1=pointc
         print "min1"
         print min1
         print "min2"
         print min2
     if min1<5:
         fet.setGeometry(QgsGeometry.fromPolyline([pts, point1]))
         fet.setAttributes(attrs)
         outFeat = QgsFeature()
         pr = layer.dataProvider()
         outGeom  = fet.geometry()
         outGeom = outGeom.combine(geom)
         outFeat.setGeometry(outGeom)
         outFeat.setAttributes(attrs)
         layer.addFeatures(([outFeat]))
         layer.deleteFeature(f.id())

     if min2<5:
         fet.setGeometry(QgsGeometry.fromPolyline([pts, point2]))
         fet.setAttributes(attrs)
         outFeat = QgsFeature()
         pr = layer.dataProvider()
         outGeom  = fet.geometry()
         outGeom = outGeom.combine(geom)
         outFeat.setGeometry(outGeom)
         outFeat.setAttributes(attrs)
         layer.addFeatures(([outFeat]))
         layer.deleteFeature(f.id())


texts:
QGIS 2.18.5

No matter what I do, if I set minimum visibility to 1,000,000 and max to 400,000 click apply and re-open properties, the max 'resets' itself to 100,000. Visa versa if I set min to 400,000 and max to 100,000,000 - the same occurs again. Is this a bug?


texts: How to zoom past 1:1433 with OpenLayers plugin in QGIS 2.6?

texts:
I am editing the polygon shapefile in QGIS 2.6 with the projection of "USA_Contiguous_Equidistant_Conic" but I am not able to zoom in the OpenLayers like Google Maps and Bing maps more than 1:1433 scale in QGIS 2.6. I verified it after installed "Tile Map Scale Plugin" as well, but still I can't able to zoom in.

Does anyone knows how to fix this?


texts: I have a list of many many rasters. I want to have a simple and easy to understand script that takes this and merges it into one big mosaic. I found a similar question but the accepted answer adivices against using R for this purpose. I could also use R but are very inclined towards Python. I really need a script as I am going to simply share it with some people that want it in this form. I am inclined to using gdal_merge.py but I wanted some thoughts on the problem first as I am potentially going to handle massive amounts of images. Why merge so many images? I also found it impractical but was something I was asked to do.


texts: Best way to find points outside state boundary that are supposed to be within specific State

texts:
I have a nationwide featureclass called bakeries. In the bakeries FC, there is a field called "State" and the field is populated with state abbreviations(AL,NJ,CT,etc). I have another feature class called United States, and that FC also has a field called "State" but it is populated with full state name.

The bakeries dataset has some poorly geocodeded points that say that they are in a state(Alabama for example) but the location of the points falls outside the state boundary, sometimes very far away.

Here is a snippet of what I have so far, you'll see that the code works but returns me everything that falls outside of the state versus the state specific data.

My goal is to identify points that say they are within a state, but in actuality, fall outside the boundary.

for row in state_rows:

    arcpy.AddMessage(str(new_name + "=" + "'" + row.getValue("STATE") + "'"))
    sqlExp = new_name + "=" + "'" + row.getValue("STATE") + "'" #making SQL Statement

    States2FL = arcpy.MakeFeatureLayer_management(States, outworkspace + "\\" + row.getValue("STATE")+ "_FL_FL",sqlExp) #Making Feature Layer of States.shp
    FC2FL = arcpy.MakeFeatureLayer_management(fc.CatalogPath,outworkspace + "\\" + row.getValue("STATE") + "_FC_FL") #Making Feature Layer of FC

    arcpy.SelectLayerByLocation_management(FC2FL,"",States2FL,"","NEW_SELECTION")#selecting all data inside respective selected state
    arcpy.SelectLayerByLocation_management(FC2FL,"","","","SWITCH_SELECTION") #switch selection

    record_countAF = int(arcpy.GetCount_management(FC2FL).getOutput(0)) #Get After Record count
    percent = (float(record_countAF)/float(record_countb4)) * 100.00 #calc percent here


texts: Dissolving only selected features in a layer?

texts: I am curious to know if it is possible to dissolve only the selected features in a layer in QGIS.  In ArcGIS this is the default behaviour (which was actually someone annoying most of the time), but it can be useful to avoid having to export all the features to be dissolved to a new file, dissolving them, and then pasting them back into the original.  Is this possible in QGIS? I do not see it in the Vector > Geoprocessing Tools > Dissolve dialog.


texts: Problems with QGIS 3D animation in Windows

texts:
When I use the 3D map view in QGIS 3.10 on Windows 10 and animate a map all the pictures turn out black. I can click the export button and get an actual picture. I have reformatted the computer and reinstalled QGIS with the same problem. It works fine on Ubuntu 18.04 with the same data. Any ideas what might be wrong.

Currently I have a Toshiba Portege M780 with Windows 10 Professional and stock video drivers.

It has been a while since I last tried to get this to work, but I think I have tried multiple computer, but all with stock drivers.


texts:
I have a filter that is currently used in URL, but it became more complicated and exceeds allowed URL length.

Is there a way to make a POST request with filter definition in a request body?

I tried to send following body, but without success. Tile is displayed, but not filtered.

<Filter xmlns:gml="http://www.opengis.net/gml">
    <PropertyIsEqualTo>
        <PropertyName>key1</PropertyName>
        <Literal>a</Literal>
    </PropertyIsEqualTo>
</Filter>


texts:
I am trying to modify the data source of a raster layer in python, then update my composer and finally save it.

To modify my data source I found what is needed here:
https://geogear.wordpress.com/2015/05/15/changing-qgis-layer-datasource-with-python/

The problem is that I lose the symbology of this layer in the legend of my composer.
I tried different codes but none update my composer legend :

iface.activeComposers()[0].composition().refreshItems()

for i in iface.activeComposers()[0].composition().items():
    if isinstance(i,QgsComposerLegend):
        legend = i
        legend.updateLegend()

The only solution I found at the moment is to use the "setAutoUpdateModel" method on my legend and switch to True / False Auto Update to force the legend's update, but this solution is not optimal.

=> I wish I could just update the symbology for a particular layer in the legend of my compose (a bit like when in the interface of my composer I click on my legend then in "property of the object" right clicks on my layer and "back to default").

Would anyone have an idea to do that?

EDIT :

I have found a solution here :
Changing legend item text with pyQGIS?

My new code :

for legendLyr in iface.mapCanvas().layers():
        renderer = legendLyr.renderer()
        myRenderer = renderer.clone()
        legendLyr.setRenderer(myRenderer)
        legendLyr.triggerRepaint()

Now my composer legend is refresh !


texts: ExportToPNG with definition query for unique values?

texts:
I have a MXD with a layer that has a field that has many values.

I need to export multiple layouts to PNG, applying a definition query to show the unique values for that layer.

i.e. Just show the features that share the same date, export, then move onto next date.

Data Driven Pages would be perfect, IF I wanted to export to PDF.

Looking for some guidance on how to script this with python or modelbuilder.  I believe a combination of ExportToPNG and applying definition queries would be the way to go?


texts: Is there a way to have a 'live' buffer in QGIS

texts:
I use QGIS (2.10.1) to create buffers around polygons.  I often create a series of buffers, and later need to edit the underlying polygons.  It is frustrating to start from scratch in re-create several buffers based on the changes when this happens.

Is there a way to create a 'live' buffer, that will automatically update to changes made to the parent polygon, or as a second best, a means to replicate with one click, the buffer settings and formatting that you originally created?


texts: How to create a polygon from a data CSV file with coordinates using QGIS?

texts:
Having the following information in a CSV file:

I want to build a polygon with all the points where the Group name is the same.

For example I want to join all the points where the group is T101 and form a polygon.

How can I achieve it using QGIS?


texts:
I am using R to develop species distribution maps.
Since I standardized my covariates during the modeling process, I needed to transform the country-wide data (i.e. riv) using the same values.

I read in my layer as a tiff with readGDAL.
I used the function attribute to obtain the following values.

However, I got the following error when I tried to transform the country-wide data.

  Error in riv - 0.19 : non-numeric argument to binary operator

Anybody know why this is not working?

riv <- readGDAL("stackkazc10a_Clip.tif") # layer to be standardized
sc.s # data set for the scaled variables
attributes (sc.s)
#$`scaled:center
riv
0.1981767
#$`scaled:scale
riv
0.2406652
riv.s <- (riv-0.1981767)/0.2406652


texts: ArcPad Checkout from python command line

texts:
I have been developing a script to automate my ArcPad checkin/checkout processes.  My problem is that I am having trouble getting arcpy to recognize the ArcPad tools from the command line.  From the example below you can see the specific error I'm getting.  There is however some recognition that the tools exist as the help window actively shows the tool and parameters, giving help suggestions while you type as it would for any other tool.  I have also followed the recommendations from this post:

Python scripting for ArcPad Checkin

I have about 100 checkin/checkouts to perform which is why I'd rather automate the process.

Example:

>>> arcpy.ArcPadCheckout_ArcPad("points_2015","","","","checkout_points.axf")
Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "C:\Program Files (x86)\ArcGIS\ArcPad10.0\DesktopTools10\ArcPad Tools.tbx",   line 66, in ArcPadCheckout
AttributeError: Object: Tool or environment <ArcPadCheckout_ArcPad> not found
>>>

Attempt from example given below, similar error message:

>>> import arcpy
>>> from arcpy import env
>>> env.workspace="H:/path_to_my_geodatabase.gdb/"
>>> arcpy.AddToolbox(r"C:\Program Files (x86)\ArcGIS\Desktop10.1\ArcToolbox\Toolboxes\ArcPad Tools.tbx")
<module 'ArcPad' (built-in)>
>>> arcpy.ArcPadCheckout_ArcPad("points_2015","","","","checkout_points.axf")
Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files (x86)\arcgis\desktop10.1\ArcToolbox\Toolboxes\ArcPad Tools.tbx", line 66, in ArcPadCheckout
AttributeError: Object: Tool or environment <ArcPadCheckout_ArcPad> not found
>>>


texts: Displaying latest satellite images and adding shapefiles on globes other than Google Earth?

texts:
Are there alternatives to Google Earth that will display the latest satellite images and can add pinpoints and are relatively intuitive?

It would be nice if I could possibly add shapefiles  on top of the satellite images.

I have tried ArcGlobe which is interesting but many layers like boundaries and lakes need to be added manually, and navigation is a little clumsy.


texts: How to set up GDAL Java Bindings for Windows with OSGeo4w?

texts:
I've been searching for ways to use java bindings for GDAL and installed the OSGeo4W packages QGIS, GDAL, and GrassGIS using the "Express Desktop Install" method.

After installing and confirming the GDAL commands functioning in the OSGeo4W Shell, I do not see the gdal-java package and all of its contents that osgeo.org says it should here. The package should include the Java SWIG bindings for GDAL/OGR and osgeo4w\lib\gdal.jar.

I have the following Windows Environment Variables set up:

 GDAL_DATA set to C:\OSGeo4W64\share\gdal
 GDAL_DRIVER_PATH set to C:\OSGeo4W64\bin\gdalplugins
 appended C:\OSGeo4W64\bin to PATH

Has the gdal-java package depreciated? (considering the web pages have been updated 7-10 years ago) If so, is there another way of installing GDAL Java bindings easily?

Computer specs: Windows 10 x64


texts:
Is it possible to retrieve a list of REST Services from an ArcGIS Server along with their MXD, and data sources (DB)??

I am using ArcGIS Enterprise 10.7.1


texts: I am trying to calculate the area of polygons in a shapefile on QGIS 2.2.  My project is in an equal area conic projection with OTF projection activated.  When I added a new field called 'AREA' and used the $area function in the field calculator, I got impossibly small numbers.  I then switched the CRS on the shapefile from the properties menu from NAD83 to the same US Equal Area Conic and tried to calculate again with the same results.  Do I need to make a copy of the shapefile and set to the CRS that way?  Or am I doing something else wrong?


texts: Hiding missing values on bar charts in QGIS

texts:
I'm trying to hide the bars when there are no values available on my bar charts. This, in order to differentiate missing values from values that are existing but are very close to zero or equal to zero. Do you know if there is a way to do this?
Here is an example on the following image, with the orange arrow indicating missing values and the blue one real values.


texts: How to delete duplicate nodes/points in QGIS?

texts:
I'm new to this site and most impressed with the help available. I have a question regarding removal of duplicate points which has been asked before.  However the QGIS solutions suggested do not appear to work for me. I'm using QGIS 2.0.1.

I have a file of about 100 attributed features which I have brought over from MapInfo into QGIS. I've run the Vector\Geometry Tools\Check geometry validity which gives about 40 duplicate points/nodes. To remove duplicates I have run the MMQGIS plugin which doesn't remove any and also the Topology Checker plugin using "must not have duplicates" but which tells me it can't find any errors.

Other suggestions have been using Generaliser plugin (which I cannot find in the repository)and also the Delete Identical (does not appear in the Data Management tools). Other suggestions are GRASS which I’m not familiar with yet.

Any suggestions from the experts as to the best/easiest way of doing this, specifically in QGIS?  Thanks

Update. The following are the steps:
1) Run Vector \ Geometry Tools\ Check geometry validity
2) Save a file with only these duplicate points
3) Select the layer with duplicates
4) Run MMQGIS\ Modify \ Delete duplicate geometries specifying the selected layer as source and a new file name for the "cleaned" layer
5) Run the geometry validity tool again on the “cleaned “ layer and get any captured duplicate points
6) Compare with original duplicates, no change

As regards an image, the only thing I get is about 40 points displaying where the duplicates are.


texts: Error Loading GeoJSON Vector in OpenLayers 3

texts:
I've followed several other posts but all the basic solutions have not worked. Basically, I can get it to show a GeoJSON correctly if I use a web link as the URL, but not a local file. For testing purposes, I will be using this GeoJSON file from OpenLayers, but the reason I need it to be a local file path is because I'm using a different GeoJSON file.

The following code works fine:

var vectorSource = new ol.source.Vector({
    format: new ol.format.GeoJSON(),
    url: 'https://openlayers.org/en/v4.6.4/examples/data/geojson/countries.geojson'
});

var vector = new ol.layer.Vector({
    source: vectorSource,
});

But, when I save that file to a local path and try linking it, I get nothing to show on the map. I do know for sure that the path exists and is correct because I have other things in the same resource directory that load fine using the same path directory:

var vectorSource = new ol.source.Vector({
    format: new ol.format.GeoJSON(),
    url: '../resources/countries.geojson'
});

var vector = new ol.layer.Vector({
    source: vectorSource,
});

What am I missing? I tried following the suggestion given here involving the IIS Manager but nothing has changed. I also tried changing the extension from .geojson to .js with no change also.


texts: OpenLayers - redrawing map after container resize

texts:
In my web application, I want to enable users to set the size of the map container.

Everything worked fine when the container was expanded slightly (apparently this is because the tiles that were just behind the border were already loaded). However, when the container is expanded significantly (in the following example, from 300 to 1000px in width), there is blank space left.

How to make the map redraw and adapt to the new size?

Calling redraw() on all layers didn't help. Neither did zooming in and out.

I tested this with the described results in Opera, Chrome and Firefox.
In IE8, surprisingly, the problem didn't occur and the map automatically adapted.

A simplified webpage for testing:

    <input type="button" id="test" value="resize">


texts: Is there a way of auto_increment for the ID column in QGIS

texts:
I'm creating a fictive map and so I need to create lots of points, lines and for sure polygons. Later on I export my data as geojson. But before that I aways have to go and give every element an unique ID.

I don't need a special sorting, like the biggest polygon gets the smallest ID or so. I just need all polygons with an ID at the end, without doing it manually like I have to do now.

Would be great if someone knows how to do that.


texts: Converting latitude/longitude into a slab number (or item number) and vice versa?

texts:
Purpose: To compress the data transmission for the longitude and the latitude. Initially these datas are in the form DD.MM.mmmm: degrees, minutes, ten thousandths of a minute. Latitude stands under 28 bits and longitude occupies 29 bits. So for each position one must transmit 57 bits.

To reduce the amount of data to be transmitted, we start from the following principle: If we set the precision of the position to 1 meter (for example) and that we divide the earth, in the considering as a sphere, in slabs of 1 square meter contiguous, we obtain a number of slabs < 2^49 so instead of transmitting latitude latitude position under 57 bits one can transmit the number of the slab under 49 bits (gain 8 bits per transmission).

Do you know how to convert longitude/latitude in a slab number (or item number) and vice vesa? what is the equation or algorithm to do it?


texts: Renaming all feature classes in file geodatabase?

texts:
I have 8 separate geodatabases in which are contained tiled feature classes all of the same name e.g. TopographicArea. So that each of my File Geodatabases contain a feature class called TopgraphicArea etc..

I would like to be able to merge the feature classes from each geodatabase of the same type (e.g all feature classes called TopgraphicArea etc..) but as it stands I can't because the names are identical.

How can I rename, or append a suffix to, each feature class so that they can subsequently be merged?

Ideally I'd like the feature classes in gdb1 to be "1" in gdb2 to be "2" and so on and preferably to be done in ModelBuilder.


texts:
I can embed a tweet in an ol3 popup using a simple iframe code stored in a geojson. The code uses a free service called twitframe.com:

<iframe border=0 frameborder=0 height=250 width=250 src=https://twitframe.com/show?url=https://twitter.com/exploreorg/status/1002995854159089664></iframe>

Working fiddle here: http://jsfiddle.net/Brian393/mr4n3w8z

However I would like to embed a full timeline, not just a single tweet from an individual's status. Instead it should be a hashtag, or even better, a moment (ie a curated collection of tweets), such as these:

https://twitter.com/hashtag/bearcam

https://twitter.com/i/moments/1006256948990312448

When I insert the URLs into the iframe above, they don't work. It seems that twitframe is designed to show only an individual tweet. Does anyone know how to embed a full timeline in an ol3 popup?


texts: ENU to rolled pitched yawed coordinate transform

texts:
I'm trying to figure out the correct way to perform a certain coordinate transformation, and I think the best way to describe it is with an analogy. Imagine a tank has rolled up the side of a hill and it must point a laser at a UFO:

The command station, knowing the UFO's GPS location as well as the tank's GPS location, informs the tank of the UFO's position via a relative East,North,Up (ENU) coordinate system:

The tank operator needs to figure out how much to rotate the turret in az and el in order to point the laser at this target, relative to its local coordinate system. I'm calling the axes Side,Forward,Roof (SFR):

In addition to having the ENU coordinates, the tank operator also knows the following information:

1) It knows the bearing of the base relative to true north (which direction it is driving when driving forward).
2) It has an inclinometer installed on the base, squared up to the F and S axis. So it gets a roll and pitch measurement. The inclinometer utilizes gravity to make this measurement, so I believe these measurements are always indicating an offset from a perfectly level plane (I believe this is the way all inclinometers work, but just thought I'd state it).

I believe all I need to do is transform the ENU coordinates to SFR coordaintes, and then do the trig to figure out the correct az and el degrees. I'm struggling a little bit on visualizing and figuring out if I should rotate for the bearing first and then the roll/pitch, or the other way around. Maybe there is some cross product needed to get the roll and pitch in terms of the ENU axis first?

Can someone help me solve this problem?


texts: Switching basemap automatically in QGIS?

texts: I want to know if there is a way to switch basemaps if some basemaps are not able to show any imagery while creating an atlas in QGIS. I understand it is possible to do it individually but when you have an atlas with hundreds of maps it would be great if this is automated.


texts: How to reproject data from [-180 180] to [0 360]?

texts:
I have both a GeoJSON and shapefile of volcanoes along the Pacific Ring of Fire. The area I want to render in the map is right around the antimeridian, so Leaflet cuts it in half. I initially tried to use VectorGrid.Slicer to wrap the data, but this was messing with the styling so I thought I would try and reproject the data itself in QGIS.

What is the best way to go about this? Duplicate the data and add 360 to each of the longitude values, or is there a more elegant solution?

I don't believe this is a repost, since most of the questions I've seen are asking how to reproject [0 360] to [-180 180].


texts: Using slope and total height for terrain analysis in QGIS?

texts:
Using QGIS I'm trying to analyse an area/DEM with these two criteria:
-Even sloping terrain steeper than 1:20 and total slope height > 5 m.

Finding areas with slope of 1:20 (5% slope) could be done with the slope function. But I can't figure out an easy way to filter these areas to only include areas with slope height over 5 meter. I.e. the vertical height form start of slope to top should be at least 5 meters.

Anyone have any ideas or thoughts on this problem?

Background: The criteria are to identity possible land slide areas (quick clay), so it's already clipped to only include areas of old marine sediments. https://www.ngu.no/en/topic/quick-clay-and-quick-clay-landslides


texts: Select By Attribute Query Using Shape Area

texts:
I have a script that assigns a rank to each feature in a feature class based on the sort order of a target field. When I use this script to assign a rank based on shape area, I cannot pull the full shape area using SearchCursor and getValue (I get 8 decimal places, instead of 13), so my queries never end up selecting any features. Is there a way to get the full float for the shape area? For now my work around is to create a new field that has the shape area rounded to 6 digits, and rank features based on that.

targetField = "Shape_Area"

#build list of unique vals
qRows = arcpy.SearchCursor(inLayer)
for row in qRows:
    uVal = row.getValue(targetField)
    if uVal not in uList:
        uList.append(uVal)
#sort list
uList.sort()

for val in uList:
    qry = '"' + targetField + '" = {}'.format(val)

# make selection based on qry
arcpy.SelectLayerByAttribute_management("inLayer_lyr", "NEW_SELECTION", qry)


texts: Automatically measure polygon length in km

texts:
I have a large number (>1000) of polygons in a shapefile. I would like to measure their approximate length. The polygons are somehow elongated, so a quantity such as the length of longest line within a polygon, or the half of its perimeter, could be good approximations for the polygons' lengths. I am new to GIS. I can load the shapefile on Google Earth Pro for example, and manually measure these lengths. But this would be a very tedious and inaccurate work. I also have access to QGIS and Python.
In short: is there a way to obtain distances (in km) between the two farthest points in a polygon automatically, or the polygon's perimeter, given a large number of polygons in a shapefile?
Preferably using Python, QGIS or Google Earth Pro.
EDIT: Since my polygons are so elongated, I would be happy with approximating the length by half of the perimeter. This could actually be more accurate than the length of longest line within, because some polygons are also curvy.


texts: Displaying time series data by changing object color in QGIS

texts:
I am looking to visualize the changes in lake chemistry over time by varying the color of the lake(s) as rendered on the map.  I have found the time manager plugin, that seems to be close to what I am looking for but I am not able to coax it to do what I want.

Ideally I would be able to put a slider on my map that would allow the user to select the date and easily move the slider back and forth to visualize the changes.

Any suggestions?


texts: Erasing extent of one layer using another layer using ArcMap?

texts:
I have two layer files which mark the extent of some vegetation types. The base layer includes some areas which does not need to be included. The second layer is solely comprised of unnecessary area coverage.

I need to remove this unnecessary area. I was hoping to subtract the unnecessary area layer from the base layer but am unsure how to do this. I don't have an advanced license so I cannot use the erase tool.

I am using ArcGIS 10.1 for Desktop.


texts: Interpreting JSON lat/lon data from an ESRI interactive map?

texts:
I got access to some data from an ESRI-based interactive map and I am trying to make sense of the data. In particular, I'd like to extract the geometry of some shapes. Below is a snippet of the JSON file. The data I'm interested in are the coordinates in the rings array.

"spatialReference": {
  "wkid": 102100,
  "latestWkid": 3857
},

...

{
  "attributes": {
    "OBJECTID": 9,
    "SEC_AREA": "Place",
    "SHAPE_Length": 0.47231263503454957,
    "SHAPE_Area": 0.0018331050756540896
  },
  "geometry": {
    "rings": [
      [
        [
          -13710176.908767894,
          6294974.636926275
        ],
        [
          -13710176.877294399,
          6294974.687023439
        ],
        [
          -13710176.612313965,
          6294974.836912473
        ],
        [
          -13710176.908767894,
          6294974.636926275
        ]
      ]
    ]
  }
}

The closest I've found is the this document which describes the rings, but does not describe the format of the coordinates themselves. They are certainly not latitudes/longitudes.

From experimentation, I've found that if I divide these coordinates by constants, one for the X and another for the Y, I do get close to the lats/lons I'm expecting, but I'd like to find out how these values are meant to be encoded. Can anyone shed some light?


texts: JTS polygon intersection calc on a geoid

texts:
I am using JTS to compute intersections between polygons with WGS84 lat/lon coordinates. I am not converting coordinates to anything 2D based, just reading in POLYGON(lat lon, lat lon) and doing geometry1.relate(geometry2).

Will that give me correct answers all the time?

My polygons are roughly 1m to 100m in size, located in US and Alaska.

If it's completely wrong, should I convert WGS84 to UTM? Using JTS/GeoTools/Proj4J ?


texts:
I am trying to align rasters in pyqgis (snap raster ArcMap equivalent).

There is a class reference in documentation, which is doing exactly what I need. QgsAlignRaster class.

Anyway, my problem is I can't get the extent right. If I put 2 rasters with different extent, like 2 overlapping rasters, the one which I want aligned will get clipped by the other, so I am missing almost half of the original raster.

Is there a way to set that I do not want  to clip it?

Documentation:
https://qgis.org/pyqgis/master/analysis/other/QgsAlignRaster.html

My code is:

def alignRaster(raster, raster_to_align, output, proj, cell_size):

 alignClass = QgsAlignRaster()
 item1 = QgsAlignRaster.Item(raster, "")
 item2 = QgsAlignRaster.Item(raster_to_align, output)

 alignClass.setRasters([item1, item2])
 alignClass.setParametersFromRaster(raster, proj)
 alignClass.setCellSize(cell_size,cell_size)
 alignClass.checkInputParameters()
 alignClass.createAndWarp(item2)

What i tried is setting the extent of raster with:

rLyr = QgsRasterLayer(raster)
clipExtent = rLyr.extent()
alignClass.setClipExtent(clipExtent)

I also tried setting clip extent to null, which is the default option I believe, since if I try to read the the clip extent, its says this:

print(alignClass.clipExtent())
>> <QgsRectangle: 0 0, 0 0>

I am out of ideas, is there a solution to my problem?

In the picture below, you can see how it gets clipped (gray raster is the new one)...


texts: Error when trying to add a new field to attribute table

texts:
I have my attribute table open, in edit mode and am trying to add a new column. the dialogue box open but I am unable to select the type" for the field and then when I press OK it gives "Failed to add field 'people' of type ''. Is the field name unique?".

I also opened the properties of the layer to try but there the add and delete columns are both grayed out.


texts: How to merge geoseries to geodataframe based on geometry attribute

texts:
I have 5868 points in a geodataframe with some columns/attributes. Between points which have distance less than 10 m, I just want to select only one point as representation in that area. I have done using the following code:

ships = gpd.read_file(r"D:\Suhendra\Riset BARATA\data ais\lego_python\kepri_201812_ship.shp")
#'ships' have 5868 data/rows. It is geodataframe with some columns

#remove the 'ships' geometry that have less than 10 m distance each other
point_nodes = list(ships['geometry'])
for i in range(len(point_nodes) - 1):
    if point_nodes[i] is None:
        continue
    for j in range(i + 1, len(point_nodes)):
        if point_nodes[j] is None:
            continue
        if point_nodes[i].distance(point_nodes[j]) < 10: #in meter
            point_nodes[j] = None

new_point_nodes = gpd.GeoSeries([node for node in point_nodes if node is not None])
#'new_point_nodes' have 5321 data, it is just geoseries with geometry information

The result is 5321 points (reduced than the original data), but it is just geoseries not geodataframe like the original data. How to do the following condition in order to get the result like the original data?


texts: Georeferencing LIDAR scans in CloudCompare?

texts:
I've got several LiDAR scans in local coordinates (point cloud with an arbitrary  0,0 as origin) that I would like to reference to real world coordinates. I've tried to Align the scan in CloudCompare using known reference locations that can be recognized in the scan, but without any success, it rotates correct, but the coordinates of the points stay in a local system.

Does anyone have experience georeferencing LiDAR scans, and can explain the workflow to be used.

Example, the real world coordinates of the five corners of the building are known:


texts: Place names on maps - why big differences between OSM, Bing and Google?

texts:
The following observations are from clicking between Chrome tabs.

Managua, Nicaragua:
At zoom 12 with Google Maps, I get about 40 place names.
At 12 zoom on OSM, I get about 17, and only two of them - Las Colinas and the airport - are the same as Google's.
At similar zoom, Bing has 11 names, and only two - the airport and Carretera a Masaya - are the same as Google's. (Yahoo's names are the same as Bing's)

So the maps are giving importance to different places.
Does anyone know the logic?

What names should be on the map of a city at zoom 12?

====

In response to the below, I'm doing maps of cities for an on-paper publication.

I always assumed that place names on maps were the major "places" in a given area, with some kind of hierarchy - like Country, state, city and City, district, suburb.

Google, OSM and Bing come up with very different place names for cities I'm interested in (leaving aside names of highways/streets).

So, I don't know how to represent these cities.
Maybe I should just take a reputable atlas and eliminate the names in small fonts.
But I'd still wonder how the mapmakers decided what should be small and what should be bigger, or there or not there.


texts: Problem with negative radius for ST_Buffer

texts:
I am very new to PostGis and I am trying to use negative radius of buffer for ST_Buffer function.

The column is geometry(MultiPolygon, 4326) imported from a shapefile. Positive radius for ST_Buffer works:

(select  ST_AsGeoJSON(geom) as no_buffer, ST_AsGeoJSON(ST_buffer(geom, 1)) as buffer from sections)

However, when I user negative radius, e.g. -1 I get:

ERROR: lwgeom_free called with unknown type (105) Invalid type
SQL state: XX000

I am sure it is some kind of noob question, but I couln't find the solution using google and StackExchange archive.


texts: What are the essential elements of a cartouche, if any...?

texts:
I'm trying to design artistic cartouches for some maps that show traditional knowledge locations from specific First Nations in Canada so the cartouche's will reflect their Nation's identity, culture, etc...

I'm not having trouble with ornamental features or colours. I struggle trying to fit more modern elements like titles, descriptions, data sources, and maybe even projection/datum information.

Are there some guidelines I should follow when trying to decide what other elements "must" be in the cartouche?

In this case, the audience will be the knowledge holder (person who contributed to the information on the map) so that they can have a hardcopy of their work to share with whom they please.

The second set of maps will be for the community showing the comprehensiveness (everyones information) of the study. Possibly displayed on a wall in the community building for example.

Below is a small image of a simple cartouche from a historical map I have for those who don't know what I mean:


texts: Creating a vector layer using the graphical modeler and a Quick OSM algorithm

texts:
I'm trying to load a vector layer via the graphic modeler into my project and the quick osm "build query inside an area" algorythm. I set up the structure like this:

Key = landuse, In = Passau, Germany.
But it does not load into my project. I am aware that i need some kind ouf output algorithm, but anything I have tried to add so far, did not work. What algorithm do I have to add after the "build query inside area" so that it acutally loads into the project?


texts: QGIS mbtiles vs xyz tiles export inconsistency

texts:
I want to use mbtiles in a project so I am generating mbtiles for a selected area of osm.

Extent (xmin, xmax, ymin, ymax):

463297.84985694464,538040.0735187809,6769150.388750425,6834675.317430764 [EPSG:3857]

I was a bit surprised to find out that I wasn't able to find a record for these values

z = 14
y = 5412
x = 8390

Query I am running:

select tile_data from tiles where zoom_level=14 and tile_column=8390 and tile_row=5412;

has no results.

However

select tile_data from tiles where zoom_level=14 and tile_column=8390;

has 28 results, where tile_row ranges from 10960 through 10987.

So I thought I would take xyz tile exports for the same area and settings to see if it has a tile for the given xyz and IT DOES!

I don't understand how this is possible! Is there something I am missing while exporting mbtiles? OR is this a bug in QGIS?


texts: I want to calculate the numbers of zoom levels that has a WMS. I had a look to the WMS spec and I guess that it has to be related with bounding box area and image size.


texts: Is it possible to convert a Geoserver SLD to Quantum QML format?

texts:
I have a web application which runs on Geoserver. I am migrating my application to desktop and am struck at applying the same color schemes and symbology. I am not able to find a way to convert SLDs to QML for implementing in desktop application.

I am aware that Quantum GIS 1.8 has the SLD input in styling. I specifically need a solution for Quantum GIS 1.6, as it is more preferred in my organization. Any suggestion is highly appreciated.


texts:
I am new to GeoNode. I installed GeoNode stable version from apt-get install geonode on Ubuntu 14, and its working properly. Now I want to change the uploaded directory to /home/user/uploaded
I changed MEDIA_ROOT path in local_settings.py and uploaded allies in etc/apache2/sites-available/geonode.conf but once a layer uploaded in base_link PostgreSQL table, the thumbs URL still stored as the original GeoNode uploaded path.
How can I change it to every new upload?
Or else is there any other way to change the uploaded directory path in GeoNode 2.4
I already searched way to change the GeoNode uploaded directory path, but I did not found a way.


texts: Python error in QGIS Ubuntu 13.10 while trying to join .csv tables

texts:
Using the MMQGIS plugin in QGIS I can't join my CSV file...

This is what QGIS returns me:

An error has occurred while executing Python code:

Traceback (most recent call last):
  File "/home/elisabet/.qgis/python/plugins/mmqgis/mmqgis_dialogs.py", line 217, in run
    joinattribute, outfilename, notfoundname, 1)
  File "/home/elisabet/.qgis/python/plugins/mmqgis/mmqgis_library.py", line 516, in mmqgis_attribute_join
    if duplication[0] != '_':
IndexError: string index out of range

Python version:
2.7.5+ (default, Sep 19 2013, 13:52:09)
[GCC 4.8.1]

QGIS version:
1.7.5-Wroclaw Wroclaw, exported

Python path: ['/usr/share/qgis/python', '/home/elisabet/.qgis/python', '/home/elisabet/.qgis/python/plugins', '/usr/share/qgis/python/plugins', '/usr/lib/python2.7', '/usr/lib/python2.7/plat-x86_64-linux-gnu', '/usr/lib/python2.7/lib-tk', '/usr/lib/python2.7/lib-old', '/usr/lib/python2.7/lib-dynload', '/usr/local/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages/PILcompat', '/usr/lib/python2.7/dist-packages/gtk-2.0', '/usr/lib/python2.7/dist-packages/ubuntu-sso-client', '/usr/lib/python2.7/dist-packages/ubuntuone-client', '/usr/lib/python2.7/dist-packages/ubuntuone-control-panel', '/usr/lib/python2.7/dist-packages/ubuntuone-storage-protocol', '/home/elisabet/.qgis/python/plugins/mmqgis/forms', '/usr/share/qgis/python/plugins/fTools/tools']


texts: Is it possible to get the bounding coordinates of a google map that a user has zoomed to using Javascript? I see that it can be done using the Google App Engine but I don't know if is possible using the Google Maps API (Version 3).


texts: Getting a list of ObjectID from a Selection Set Object using ArcObjects 9.3 in VB.Net

texts:
I am trying to get back a list of ObjectIDs from a selection set using VB.Net in ArcObjects for ArcGIS Desktop 9.3.  I have been passing an object populated using the ".IDs" method  from the ISelectionSet interface ((which returns an IEnum2 object with the OBJECTIDS for the selected features) to the following function:

    Public Function iterateOIDS(ByRef OID_Enumerator As ESRI.ArcGIS.Geodatabase.IEnumIDs)
    Dim itemID As  = OID_Enumerator.Next()
    Dim listofOIDS As List(Of Integer)

    Do While itemID <> -1
        listofOIDS.Add(itemID)
        itemID = OID_Enumerator.Next()

    Loop

    Return listofOIDS

End Function

However, this crashes ArcGIS and results in Visual Studio throwing an LoaderLock error, if its in Debug Mode.  The offending piece of code is apparently the "listofOIDs.Add(itemID)" inside the While Loop.

What am I doing wrong?  Any light that anyone can shed on this problem would be awesome!


texts: Issue when reprojecting shapefile using OSGEO in Python: prj file missing and kernel restart

texts:
I tried to reproject a shapefile based on the projection of a raster file following this example (section "Reproject a Layer").
My code is as follows:
    # open the raster tif and extract crs
    t = "path/to/raster.tif"
    tif = gdal.Open(t)
    crs1 = tif.GetSpatialRef()#.GetProjection()
    tif = None
    # open shapefile
    drv = ogr.GetDriverByName("ESRI Shapefile")
    s = "path/to/shapefile.shp"
    # set name of the new shapefile
    sUTM = "shapefile_UTM.shp"
    LayerName = "shapefile_UTM"
    if os.path.exists(sUTM):
        drv.DeleteDataSource(sUTM)
    shp = drv.Open(s)
    layer = shp.GetLayer()
    crs0 = layer.GetSpatialRef()
    # transformation from crs0 to crs1
    coordTrans = osr.CoordinateTransformation(crs0, crs1)
    # create output shapefile
    outSHP = drv.CreateDataSource(sUTM)
    outLayer = outSHP.CreateLayer(LayerName,
                                  geom_type = ogr.wkbMultiPolygon)
    lyrdef = layer.GetLayerDefn()
    for i in range(0, lyrdef.GetFieldCount()):
        fieldDefn = lyrdef.GetFieldDefn(i)
        outLayer.CreateField(fieldDefn)
    outdef = outLayer.GetLayerDefn()
    feature = layer.GetNextFeature()
    while feature:
        geom = feature.GetGeometryRef()
        geom.Transform(coordTrans)
        outFeature = ogr.Feature(outdef)
        outFeature.SetGeometry(geom)
        for i in range(0, outdef.GetFieldCount()):
            outFeature.SetField(outdef.GetFieldDefn(i).GetNameRef(),
                                feature.GetField(i))
        outLayer.CreateFeature(outFeature)
        outFeature = None
        feature = layer.GetNextFeature()
    layer = None
    shp = None
    outSHP = None

However, when I run this code, it seems to run into some issues, because it ends with Restarting kernel...  and all variables are gone.
Moreover, the output shapefile has only .dbf, .shp and .shx files. There is no .prj file. Of course, I could write this file in an additional step, but I think there is some underlying issue somewhere which I would prefer to resolve instead of fixing the outcome.
What is the error in the above shown code?
Edit:
The reference system I want to transform the shapefile to—i.e., what print(crs1) gives—is:
PROJCS["WGS 84 / UTM zone 35S",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",27],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32735"]]


texts: Location data for Nuclear Power stations in France

texts:
As it says on the title, does anyone know a place I can look at and/or download location data for Nuclear power installations in France?  I believe it is well known data, but I cannot find it anywhere.

It's a very innocuous request, based on choosing where to buy a new pad in France and is nothing sinister.

Thanks in advance.


texts: Getting clip raster ERROR 1: Cutline feature without geometry?

texts:
I am using QGIS 2.18.15. I tried clipping a raster

The mask layer is a shapefile. And it comes out as "ERROR 1: Cutline feature without a geometry"

I checked the point data that I interpolated and it's in the proper coordinate referencing system.

The intention is to clip the raster to the map at the back.


texts: TileStache Dynamic query based on zoom level

texts:
I wanted to make tilestache to be able have a different query for different zoom level.
For certain zoom level I wanted TileStache to query from different postgres table.
How can I achieve this? Is this even possible in TileStache?


texts: Enable vector layer using an expression via PyQGIS

texts:
How do I go about enabling a vector layer using an expression via PyQGIS?

See screenshot below for location of field, and a sample expression:

I'm new to using PyQGIS and not sure if I should be using QgsSymbolRenderContext, but I can't figure out how to get it working.

Here's what I've tried without success:

qrc = QgsRenderContext()
u = QgsUnitTypes.RenderUnit()
expression = QgsExpressionContextScope("@atlas_pagename = 'Column_2'")
symbolScope = QgsSymbolRenderContext(qrc,u)
symbolScope.setExpressionContextScope(expression)


texts: ArcGIS for INSPIRE, build Protected sites Inspire layer

texts:
Can anyone add some infomation about building protected sites inspire layer?
I cannot understand the whole specification schema and I didn't know how to use it in a real database.


texts: Splitting road network into non overlapping segments in QGIS

texts:
I am having a road network shapefile. I want to split it into smaller edges in such a way that each edge consist of two diversions. I am unable to do so.

Here n1, n2 and n3 are any three adjacent points in the road network which are diversions.

I want to convert this road network into a graph.

Image 1: n1-----n2 (desirable output)

Image 2: n1----n2----n3 (overflow condition/ undesirable)

Image 3: n1-----*******n2 (underflow condition/ undesirable)


texts:
I need to create a series of points along a line at a specific interval. I am using ArcMap 10.2.1 and the obvious tool appears to be Construct Points.

I have tried to follow the instruction steps given on the ESRI website. However, the Edit tool on the editor toolbar is not selectable so have used Editor > Start Editing and selecting the line feature instead.

At no time is Construct Points selectable on the Editor drop-down menu even after trying various thing. What am I doing wrong?

Apologies for the potential simplicity of this problem but I just can't find a solution!


texts:
Does anyone know of scripts for merging TIGER/Line shp files? It is for extracting the road features only, but for the whole of USA.

They are provided as shp files for one county at a time here:
http://www.census.gov/geo/www/tiger/

If anyone sells the data, a link is also appreciated.


texts: Python QGIS : Adding features to a temp vector layer with their attributes

texts:
I try to create a temporary vector layer with data from another vector layer that are inside a polygon.
My code works perfectly except that in the temp vector layer, the fields of the attribute table are all empty and I can't figure how to fix this.
The fields of my temporary vectory layer are the same as in the layer with my data.

def site(site):
 import processing
 #Creating the temp vector layer
 df_site = QgsVectorLayer('Polygon?crs=epsg:2154', 'polygon' , 'memory')
 #Selection of the polygon
 perim = QgsProject.instance().mapLayersByName('perimetre_exposition')[0]
 it = perim.getFeatures( QgsFeatureRequest( QgsExpression( "\"Site\"=' {}'".format(site))))
 ids = [i.id() for i in it]
 site = perim.getFeature(ids[0])
 perim_site = QgsVectorLayer('Polygon?crs=epsg:2154', 'polygon' , 'memory')
 perim_site.dataProvider().addFeature(site)

 df_site.startEditing()
 #Vector layer including my data (named bats)
 bats = QgsProject.instance().mapLayersByName('batiments_landuse_48sites')[0]
 #Adding the fields from the layer with my data to the temp layer
 for i in bats.dataProvider().fields():
     df_site.addAttribute(i)

 bats.removeSelection()
 processing.run("native:selectbylocation",{'INPUT': bats,'PREDICATE': 0, 'INTERSECT': perim_site, 'METHOD':0})
 bats_site = bats.getSelectedFeatures()
 df_site.dataProvider().addFeatures(bats_site)
 del(bats,bats_site)
 QgsProject.instance().addMapLayer(df_site)


texts: Extract domain from ArcGIS Online feature service using Python API

texts:
I have a feature service with several coded value domains that I'm accessing as a spatially enabled dataframe using the python API.

As far as I can tell, the dataframe only holds the code for each field. I'd like to  update those to the description from the domain. Is there a way to access the domain of the feature service through the API?


texts: Changing active layer in combo box in PyQgis plugin?

texts:
I came across with a problem. I don't know how to change active layer by using a combo box object (QComboBox). I manage to create a list of names but I don't know how to change active layer in QGIS plugin Combo Box yet.

index = QgsMapLayerRegistry.instance().mapLayers().values()

for i in index:
    if i.type() == QgsMapLayer.VectorLayer:
        self.dockwidget.active.addItem( i.name(), i )

Above code creates a list of my layers, but I can not change active layer with this method. Whenever active layer is changed, combo box should also update its reference.

Can you help me? I am using Qgis 2.16.


texts: How can I export MBTiles from TileMill 2?

texts:
In TileMill 1 it was easy to export MBTiles for a selected region with my own style. In TileMill 2 I can only package the style or export vector tiles from a previously imported datasource.

How can I export MBTiles with my custom style? Thanks for any hints. I haven’t been able to Google anything relevant.


texts:
This code will go to WATERLINES.gdb, go through the gdb to see if the fc exist (based on a wildcard query), then add the fc to the mxd as desired.  My questions is if the fc doesn't exist, the "else" isn't met.  That is, the warning "No Waterline data exist for the PWS." isn't added to the tool dialog.  It just keeps running the rest of the code.  I'm curious as to why?

Edit: I added more of the code.

# Import modules
import arcpy
import os

# Set overwrite option
from arcpy import env

env.overwriteOutput = True

pws = arcpy.GetParameterAsText(0)

waterlines = "WATERLINE_" + pws + "_*"
wl_layer = "Waterlines " + pws
wl_symb = r"G:\Source Water Protection\PWS_ID_Tool\WATERLINES.lyr"

# Set current mxd and dataframe.
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "")[0]

env.workspace = r"G:\WATER_FACILITY\WATERLINES.gdb"
wlFclass = arcpy.ListFeatureClasses(waterlines)
arcpy.AddMessage(wlFclass)

for wl in wlFclass:
    wlPath = os.path.join(env.workspace,wl)
    if arcpy.Exists(wlPath):
        arcpy.AddMessage("Waterline data exist for this PWS.")
        arcpy.MakeFeatureLayer_management(wlPath, wl_layer)
        arcpy.ApplySymbologyFromLayer_management(wl_layer, wl_symb)
        newLayer = arcpy.mapping.Layer(wl_layer)
        arcpy.mapping.AddLayer(df, newLayer, "TOP")
    else:
        arcpy.AddWarning("No waterline data exist for this PWS.")

Edit: the wildcard creates a list of one fc and will only ever have one. to go from unicode to a "regular" string. then create the path for the one fc to use in the exist function. when I run the tool when a fc exist i get this messgae in the tool dialog:

Executing: waterlinestest 465
Start Time: Thu Sep 10 12:33:57 2015
Running script waterlinestest...
[u'WATERLINE_465_CAW']
Waterline data exist for this PWS.
Completed script waterlinestest...
Succeeded at Thu Sep 10 12:34:45 2015 (Elapsed Time: 47.36 seconds)

if there is no fc:
Executing: waterlinestest 492
Start Time: Thu Sep 10 12:41:26 2015
Running script waterlinestest...
[]
Completed script waterlinestest...
Succeeded at Thu Sep 10 12:42:13 2015 (Elapsed Time: 46.85 seconds)

So yes, even though the list created has nothing in it, it recognizes that it "exists".


texts: Writing PyQGIS script to modify layer (add field and fill it) inside QGIS?

texts:
I have the following PyQGIS script with QGIS 3 that tries to ask for a layer, add a new field in that layer and modify the value of this field. But the code doesn't work. I get either Seems there is no valid script in the file when I use layer = QgsVectorLayer(Layer) or AttributeError: 'NoneType' object has no attribute 'isValid' when I use QgsProcessingUtils.mapLayerFromString :

##Layer=vector
context = QgsProcessingContext()
layer= QgsProcessingUtils.mapLayerFromString(Layer, context)

layer.startEditing()

if layer.dataProvider().fieldNameIndex("new_col_name") == -1:
    layer.dataProvider().addAttributes([QgsField("new_col_name", QVariant.String)])
layer.updateFields()

id_new_col= feuilles.dataProvider().fieldNameIndex("new_col_name")

for feature in layer.getFeatures():
  column_value= feature["column"]
  ... some operations ...
  layer.changeAttributeValue(feature.id(), id_new_col, new_column_value)

layer.commitChanges()


texts: Is there any Map event fired when map is loaded completely when a device is oriented from portrait to landscape mode?

texts:
We are using ESRI ArcGIS JavaScript API in one of our Mobile web applications. When the map is loaded for the first time an Onload event is fired. But when we change orientation from Portrait to Landscape, the OnLoad event does not get fired. We need such event so that when a map is completely loaded in Landscape mode we can perform some operation.

We know that there is Orientation change event in JavaScript but we need an event which should fire when the map is loaded completely in Landscape mode after orientation change event is fired.


texts: Why my 'add basemap' data seems taking too long to upload the maps in ArcGIS?

texts:
Can i know the easiest way or any other way to upload Big maps hybrid from add basemap in ARcGIS 10? because it seems taking quite a long time to upload the map like around more than 1 hour and still nothing happen. The internet connections is excellent and the desktop use are high speed. Is it because it still no comportable with something? can i know the best other way to add basemap from arcgis into layer than add data from GE to match the coordinate with other layer??

Thank you


texts: Where should I add the API key in ORS Tools?

texts:
I want to make an isochrone for my halte and station points. But I have to configure ORS Tools with my API key. The problem is, I can't find any configuration button in my ORS Tools plugin.
How should I add my API key to my ORS Tools?
My QGIS view:


texts: Exporting georeferenced image (jpeg, png) to .dxf file in QGIS?

texts:
Is it possible to export a georeferenced image, such as a jpeg with a world file, to a .dxf file using QGIS version 2.14.21?

The resulting .dxf file would only be the georeferenced image with no underlying information to be used as a basemap in another program. The dxf export option only allows vectors to be exported and I have not seen any plugins that would help. My company used a program that was able to do this but we no longer have access to it so I was wondering if this was possible in QGIS.


texts: Iterate cost distance/cost path over multiple points

texts:
I have multiple source points in one shapefile and one destination point in a separate file. I would like to find the least cost path from each source point to the destination point.  My costs are weighted using a DEM (elevation only).  I am using ModelBuilder to iterate over each source point to produce a cost distance and cost backlink raster to feed into the cost path tool.

I've run the iterator and cost distance tool and they function properly, but when I add in the cost path tool, I get multiple errors after the first source point run:

ERROR 010236: Error in building VAT.
ERROR 010067: Error in executing grid expression.
ERROR 010414: Error in creating a Raster from a geodataset.
Failed to execute (Cost Path).

I'm not very Python-savvy and am using ModelBuilder.


texts: Using external data source with Leaflet / GeoJSON?

texts:
First of all, I've been looking everywhere and I have not found a functional solution.

I need to display the points from a data file exported from Qgis.

I tried Qgis2web, but I'd like to put everything up in the Leaflet from the beginning and then just update the data from Qgis.

Explorer does not report any errors, but it does not show points with markers.

Here is index:

 <!doctype html>

      <meta charset="utf-8">
        <meta http-equiv="X-UA-Compatible" content="IE=edge">
        <meta name="viewport" content="initial-scale=1,user-scalable=no,maximum-scale=1,width=device-width">
        <meta name="mobile-web-app-capable" content="yes">
        <meta name="apple-mobile-web-app-capable" content="yes">
        <title>PinLeX</title>

        <link rel="stylesheet" href="css/leaflet.css" />

And here map set (leafletmap.js):

var map;

  function initmap() {
    // set up the map
    map = new L.Map('map');

    // create the tile layer with correct attribution
    var osmUrl='http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png';
    var osmAttrib='Map data © OpenStreetMap contributors';
    var osm = new L.TileLayer(osmUrl, {
  minZoom: 8,
  maxZoom: 18,
  attribution: osmAttrib
  })

    // start the map in defined location
    map.setView(new L.LatLng(48.974244, 14.478729),9);
    map.addLayer(osm);

  // PinLayer
  var pinpoint = L.icon({
  iconUrl: 'markers/pinpoint.svg',
  iconSize: [26.599999999999998, 26.599999999999998],
  iconAnchor: [16, 37],
  popupAnchor: [0, -28]
  });

  L.geoJson(PinLayer, {
  pointToLayer: function (feature, latlng) {
            return L.marker(latlng, {icon: pinpoint});
    },
  }).addTo(map);
}

And here GeoJson (transformed to js - PinLayer.js)

var PinLayer = {
"type": "FeatureCollection",
"name": "PinLayer",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::3857" } },
"features": [
{ "type": "Feature", "properties": { "id": null, "name": "Senovážné náměstí", "country": "Česká republika", "desc": "Popis není k dispozici.", "status": "v provozu", "photocell": "ne", "admin": "Statutární město České Budějovice (?)", "update": "7\/4\/2018", "spec": "České Budějovice", "type": "standartní pítko", "photo": null }, "geometry": { "type": "Point", "coordinates": [ 1611753.67562, 6270485.80685 ] } },
{ "type": "Feature", "properties": { "id": null, "name": "Katedrála sv. Mikuláše", "country": "Česká republika", "desc": "Popis není k dispozici.", "status": "v provozu", "photocell": "ne", "admin": "Statutární město České Budějovice (?)", "update": "12\/4\/2018", "spec": "České Budějovice", "type": "standartní pítko", "photo": null }, "geometry": { "type": "Point", "coordinates": [ 1611428.44996, 6270639.98718 ] } }
]
}


texts: Custom ArcGIS Server Web Service can't find Esri assemblies

texts:
I have ArcGIS Server 10 installed on a server, with Web ADF Runtime and Web Applications installed, and I have a custom web service that uses a couple of ESRI.ArcGIS assemblies deployed on that server. When I run the web service, I get the following error:

  Could not load file or assembly 'ESRI.ArcGIS.Geodatabase,
  Version=10.0.0.0, Culture=neutral, PublicKeyToken=......" or one of
  its dependencies. The system cannot find the file specified.

The web service uses .NET 3.5 and the ESRI assemblies are registered in the GAC. I also tried including them in a Bin folder in the app with no luck. Any help is much appreciated.


texts: Exporting Large Number of Landsat Images from Google Earth Engine

texts:
I am interested in downloading tiff images of cloud-free landsat data from around 100,000 areas of interest in Google Earth Engine. I have written the export script but the process of exporting takes a long time. I am doing this on the web code-editor. I want to know if there is a faster way of doing this. Will offline downloading using python ee package work faster?

I am interested in downloading from GEE as it contains processed clean data with clips and other such functionality and hence am not looking to use directly available data perhaps on USGS site which would also require me to use Qgis and do much of the filtering and pre-processing on my own...


texts: Saving raster as *.tif with tfw worldfile

texts:
I am using QGIS 2.16 an I need to save a raster file as *.tif with its world file as *.tfw but I only manage get *.aux.xml world files? Any ideas?

Ok, I just found out that when i "export as image" it creates a *.tfw. Minor problem seems to be that i am restricted to the screen resolution.


texts: Arcgis, Convert analyse result in numeric value

texts:
I have a problem that seems really simple but blocking me... in Arcgis i have script python where i want to extract the resolution of a raster and use this result like numeric value for use in a mathematical function.

But when i try to use it, here i have my problem because it is not a numeric value. So i want to convert the result of my request (Get raster Properties) in numeric value. Some body have a solution?

      a = arcpy.GetRasterProperties_management(r"Topographie\Raster d'élévation \Lidar\Elevation LastReturns.tif","CELLSIZEX")
      Result: Cellsize in x direction = 1.000000

      type(a)
      Result: class 'arcpy.arcobjects.arcobjects.Result'

and if i try :
          File "<string>", line 1, in <module>
          TypeError: float() argument must be a string or a number


texts: Manipulating polygon data imported from Google BigQuery GIS within R

texts:
I have some polygon shapefiles (representing UK postcodes) which I have uploaded to BigQuery, using this guide, hoping to take advantage of the inbuilt GIS support described here.

I am trying to import these into R, and from there plot these in an R Shiny app (or use for other analysis).

Whilst I have managed to import the polygon data into R, using the "bigrquery" package, the data extracted is in an unfamiliar format to me - a table where the first column is a GeoJSON string for each shape (containing only the geography), and the other columns are attributes:

# A tibble: 5 x 4
  geom                                                            POSTCODE   UPP PC_AREA
  <chr>                                                           <chr>    <int> <chr>
1 "{\"type\":\"Polygon\",\"coordinates\":[[[-2.095015148797922,5… AB10 1AL     3 AB
2 "{\"type\":\"Polygon\",\"coordinates\":[[[-2.095031775835887,5… AB10 1AN     4 AB
3 "{\"type\":\"Polygon\",\"coordinates\":[[[-2.096682747887519,5… AB10 1AP     5 AB
4 "{\"type\":\"Polygon\",\"coordinates\":[[[-2.09786728657654,57… AB10 1AS     6 AB
5 "{\"type\":\"Polygon\",\"coordinates\":[[[-2.095637704208373,5… AB10 1AU     7 AB

I am trying to convert this within R to something more familiar e.g. an sf object or a SpatialPolygonsDataFrame, which will be easier to further manipulate/plot etc.

Is there a function I can use to do this?


texts: Detecting trends in local business web services

texts:
Do you know any research paper about spatial data mining algorithms to retrieve trends about local business services (such as google places.. etc)

let's say most rated restaurants this week.. etc ?

I need scientific material about it.


texts: Can I create tiles in Qgis and upload them to Arcgis Online

texts: I have a large shapefile that I want to convert to tiles. Now I need to load those tiles up to Arcgis Online, but I only have Qgis 2.8 desktop--since I work on Linux. I know that Arcgis desktop can load tiles directly to Arcgis online. But, am I able to upload tiles created on Qgis to Arcgis Online?


texts: Autoregressive parameter from autocorrelation using R

texts:
I was wondering how to get autoregressive (AR) parameter ("rho") from an autocorrelation plot using R. Below I generated a correlated data and then provided a simple approach to produce an autocorrelation plot for the data.

#load packages
require(spatial)
require(MASS)
library(gstat)
require(rgdal)
require(sp)
#grid
xy <- expand.grid(1:100, 1:100)
names(xy) <- c("x","y")

# create spatially autocorrelated data
g <- gstat(formula=z~1, locations=~x+y, dummy=T, beta=1,
           model=vgm(psill=0.1, model="Exp", range=10), nmax=10)

# simulate
yy <- predict(g, newdata=xy, nsim=1)
coordinates(yy)=~x+y
spplot(yy, col.regions = rainbow(7)) #plot

#derive autocorrelation structure
a<-data.frame(yy)
colnames(a)<-c("x","y","z")
a <- surf.ls(2, a) #fits a trend surface by least-squares
cc<-correlogram(a, 100, xlim=c(0,100), ylab="Autocorr",
                xlab="Distance") #plot correlogram
cc #details

The last plot is indicative of an AR-model- the plot indicates autocorrelation coefficient as a function of distance, d; what I'm seeking is the autoregressive parameter ("rho" (values between -1 and +1)) from this AR model, which is indicative of the strength of autocorrelation. The relation between "rho" and correlation coefficient is defined by the following (self-explanatory) equation:


texts: Constructing a graph from a city in Python using Overpass

texts:
I'm looking to construct something exactly like this person
Building a graph from OSM?

That is I need to construct a graph with intersections as nodes and roads as edges for a simulation. Unlike this user I am doing this in Python and not PostgreSQL.

I assume I am going to need to use the Overpass Python API to get this data but I have never worked with this API before. I am competent enough in Python that as long as I can figure out how to request this data I will be fine.

However, I cant find any intuitive examples on how to retrieve this data in Python.