texts: Upper and lower bounds from gstat::krige() and gstat::idw()

texts:
Is it possible to obtain lower and upper bound values from a kriging estimation using gstat::krige() in the same way it can be obtained from a simple lm() model?
The following lm() example comes from http://pebesma.staff.ifgi.de/modellierung/course.pdf (end of page 22, 23)
library(sp)
library(gstat)
data(meuse)
names(meuse)
data(meuse.grid)
names(meuse.grid)

#Now try:
coordinates(meuse) = c("x", "y")
coordinates(meuse.grid) = c("x", "y")
gridded(meuse.grid) = TRUE

# We can use the prediction errors provided
# by lm to calculate regression prediction intervals:
plot(log(zinc) ~ sqrt(dist), meuse)
x = 0:100/100
lzn.lm = lm(log(zinc) ~ sqrt(dist), meuse)
pr = predict(lzn.lm, data.frame(dist = x), interval = "prediction")
abline(lzn.lm)
lines(sqrt(x), pr[, 2], lty = "dashed", col = "red")
lines(sqrt(x), pr[, 3], lty = "dashed", col = "red")
pr = predict(lzn.lm, data.frame(dist = x), interval = "confidence")
lines(sqrt(x), pr[, 2], lty = "dashed", col = "blue")
lines(sqrt(x), pr[, 3], lty = "dashed", col = "blue")
title("95% confidence intervals (blue) and prediction intervals (red)")

#Next, we can make maps of the 95% prediction intervals:
pr = predict(lzn.lm, as.data.frame(meuse.grid), interval = "prediction")
meuse.grid$fit = pr[, 1]
meuse.grid$upper = pr[, 3]
meuse.grid$lower = pr[, 2]
spplot(meuse.grid, c("fit", "lower", "upper"), sp.layout = meuse.grid, layout = c(3, 1))

I want to obtain these fit/lwr/upr values inside my kriging output object, next to 'var1.pred' (i.e., kriging_output$lwr). Is it possible? And can the same be done using gstat::idw()?
I tried kriging_output = predict(gstat_model, new_data, interval = "prediction") but it only returns 'var1.pred', not fit/lwr/upr.
Not sure if it's relevant but my gstat_model involves four predictors (DA ~ temp+Temp.Seasonality+um.rel+CWD)


texts: How can I hide a tool from command window under category

texts:
I have created a tool, a button and a dockable window under one project. To open/close the dockable window, button is used. For some specific work, by clicking some specific button from dockable window, I have activated the tool. So, basically the tool is used from the dockable window. So, the tool shouldn't visible to the user in the command window under the category.

I couldn't find any document to hide the tool. Is there any process that can help me to hide the tool?


texts: Where can I find official Land Use/Land Cover GIS data for Ecuador?

texts:
Right now I'm specifically looking for Land Use/Land Cover data but things like roads and soils would be handy in the future as well. In the U.S. my source for nationwide LULC data would be the Multi-Resolution Land Characteristics (MRLC) Consortium but I don't know the equivalent group/agency in Ecuador, if there is one.

I do have some LULC data for Ecuador and it is pretty detailed, so I think that there is official LULC data for the country but I need to be able to cite my sources. Also, the data I have (from 2008) is in vector format but it was obviously converted from raster, so it would be great if I could get a hold of the original, untouched raster and its metadata. Basically, I'm trying to identify the provenance of the data I do have, and find more of it (for different years) if it exists.


texts: Bivariate Mapping (running a single variable regression and mapping r-squared value)

texts:
I am hoping to produce a bivariate map displaying the correlation between two variables (specifically percent SNAP participation by county against several social variables). I would like to produce something similar to the bivariate maps produced in Geoviz Tool Kit, that can then be exported in a finished format with a scatter plot and map of R-squared values.

I have looked into GeoDa, but have been unable to figure this out.

Attached are examples of what I hope to produce from Geoviz ToolKit.


texts: Preserving Field Name & Alias in Select() function

texts:
I am unable to get an output using arcpy.analysis.Select() where the field aliases and names are preserved.

from arcpy.analysis import *

Select("Input.shp", "Output.shp", "FLOOR= 0")

The codeblock above gives columns names based on the table name such as Input_1, Input_2, Input_3, etc... is there a workaround for this?


texts:
I have this information for a custom CRS.

PROJCS["NAD_1983_HARN_Adj_MN_Hubbard_Feet",GEOGCS["GCS_North_American_1983_HARN",DATUM["D_North_American_1983_HARN",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500035.7771574167],PARAMETER["False_Northing",100007.154316],PARAMETER["Central_Meridian",-94.92055555555557],PARAMETER["Scale_Factor",1.000071553661],PARAMETER["Latitude_Of_Origin",46.80361111111111],UNIT["Foot_US",0.3048006096012192]]

How would I enter it in urs code?


texts: Point intersect-setting tolerances in ArcGIS Desktop?

texts:
I need to find out which polygon a series of points fall in.  Normally I would do this with an intersect so I can then export the data and process it externally.  As the point layer contains some coastal features some points fall outside of the polygon boundaries.

Is it possible to setup a tolerance such that if a value is going to return as null (or falls within a certain distance of a boundary) it returns the value of the nearest polygon?


texts: Standalone QGIS application - handling QgsVectorLayer load errors

texts:
I am creating an application with QGIS libraries, and in the event that the app can't add a specified layer (e.g. no connection to host, no password specified etc), it causes the app to crash, although I do get some output in the console, like QGIS displays in a message box:

credentials for dbname='testdb' host=testserver port=5432
message: fe_sendauth: no password supplied

This seems to be returned at the point where I create my QgsVectorLayer instance because the app becomes unresponsive and doesn't progress further.

My script to load the layer:

def add_postgis_vector_layer(self, username, password, layer_name, label, geometry, id):
    uri = core.QgsDataSourceURI()
    uri.setConnection(HOST, PORT, DATABASE, username, password)
    uri.setDataSource(SCHEMA, layer_name, geometry, '', id)
    layer = core.QgsVectorLayer(uri.uri(), label, 'postgres', False)  # Error message occurs at this point
    if layer.isValid():
        self.canvas.setExtent(layer.extent())
        self.canvas.zoomScale(self.canvas.scale() * 1.1)
        self.canvas.setWheelAction(2)
        # Set canvas map layers
        layers = [gui.QgsMapCanvasLayer(layer)]
        self.canvas.setLayerSet(layers)

The QApplication subclass:

class myApplication(QtGui.QApplication):
    def __init__(self, args):
        super(myApplication, self).__init__(args)
        core.QgsApplication.setPrefixPath(qgis_prefix, True)
        core.QgsApplication.initQgis()

How do I access this connection error? I can't call layer.error() (not sure if this would even work) because the app has crashed already and won't process any more code.

Edit:
I'm thinking I could avoid the error by performing a query using psycopg2 with the parameters I am about to use. Still doesn't address the issue of retrieving any errors in the connection.


texts: How to combine data from overlapping polygons?

texts:
Here's what I'm working with:

Shapefile of the current 32 Congressional districts for all states
Shapefile of the tract boundaries for Texas
A DBF format table of Hispanic counts and proportions for Texas census tracts (joined to the tract boundaries)

Using this data set, I'm trying to determine the following based upon the centers of the census tracts:

Total Hispanic population for each Texas Congressional district
Hispanic proportion of total population for each Texas Congressional district
Count of majority Hispanic census tracts for each Texas Congressional district
The shape index value for Texas Congressional district

I'm also trying to create a map applying shading symbology to display the number of majority Hispanic census tracts for each Texas Congressional district.

My problem arises when I try to combine the attribute data from the census tracts with the congressional districts. I need to combine all of the tracts within individual congressional districts and sum up their Hispanic populations. Dissolve isn't working, because the congressional districts and census tracts don’t share any common attributes. When I join the shape files spatially, most of the attribute data gets lost.


texts:
None of the threads in stack exchange related to this error helped me to resolve this issue. Hence I am posting it.

I have installed postgres 9.3 and PostGIS 2.1
Now I have done the following

su postgres
createuser -U postgres gis
createdb -U postgres -E UTF8 -O gis gis
createlang -U postgres plpgsql gis

psql -U postgres -d gis -f /usr/pgsql-9.3/share/contrib/postgis-2.1/postgis.sql
psql -U postgres -d gis -f /usr/pgsql-9.3/share/contrib/postgis-2.1/spatial_ref_sys.sql

Now when I try to create extension for postgis with the following

psql --username=postgres --dbname=gis -c "CREATE EXTENSION postgis;"

I get this:

ERROR:  PostGIS is already installed in schema 'public', uninstall it first

gis=# \dt
              List of relations
 Schema |      Name       | Type  |  Owner
--------+-----------------+-------+----------
 public | spatial_ref_sys | table | postgres
(1 row)

gis=#

This is the table

I tried altering the db with search path as well but no luck  ALTER DATABASE gis SET search_path = postgis;

This is the list of my dbs

gis=# \l
                                  List of databases
   Name    |  Owner   | Encoding |   Collate   |    Ctype    |   Access privileges
-----------+----------+----------+-------------+-------------+-----------------------
 gis       | gis      | UTF8     | en_US.UTF-8 | en_US.UTF-8 |
 postgres  | postgres | UTF8     | en_US.UTF-8 | en_US.UTF-8 |
 template0 | postgres | UTF8     | en_US.UTF-8 | en_US.UTF-8 | =c/postgres          +
           |          |          |             |             | postgres=CTc/postgres
 template1 | postgres | UTF8     | en_US.UTF-8 | en_US.UTF-8 | =c/postgres          +
           |          |          |             |             | postgres=CTc/postgres
(4 rows)

gis=#


texts: Is it possible to align Text or other drawn graphics to a feature class?

texts:
I regularly have to input measurements such as 15' or drawn lines (not feature class lines, just a graphic) to building footprints that are polygon FC's.

Are there any tools to line up/align graphics to feature classes in Data View?


texts:
When reading a las file I'm getting the following error:

> readLAS("/tmp/chunk89/u_6300085100_2015.las")
Error: Invalid data: ScanAngle greater than -180 degrees.
In addition: Warning messages:
1: There are 1091 points flagged 'withheld'.
2: There are 1091 points flagged 'synthetic'.

Is there anything that can be done to fix this?


texts: Hiding invisible map layer in legend using QGIS

texts:
I am working with QGIS 3.14 and 3.16 on Mac and Windows and I am really struggling in finding a way to remove disabled map layers in the legend view of the Layout Composer. In prior versions there was this "Filter legend by map content" button (see here: https://docs.qgis.org/3.4/en/docs/user_manual/print_composer/composer_items/composer_legend.html#legend-items) but this option seems to be removed for any reason (see here https://docs.qgis.org/3.16/en/docs/user_manual/print_composer/composer_items/composer_legend.html#legend-items)
The "Filter legend by map content" button is still available from the layer panel but this is not working for the legend in the layout. I also tried the legend item checkbox "Only show items inside linked map" without any success.
Are there any ideas or hints how to solve this issue?
Maybe it can be filtered using an expression - but how is the script looking that filters all visible/enabled map layers for the layout legend
I am wondering because usually you have a lot of work data inside a GIS project and want to display only few of them on a map, so I think I am not alone with this issue.


texts: Install both postgis and postgresql-9.1-postgis packages?

texts:
I've installed the PostgreSQL package (ver. 9.1) on a Debian 7 server using the command "apt-get install postgresql postgresql-contrib".  I now want to install PostGIS as well.  Should I install both the "postgis" package and the "postgresql-9.1-postgis" package or just the latter?  It's not clear from the descriptions on their Debian wiki pages what the differences are between the two packages.

I ask because I've seen some blog posts where the author has installed both packages but yet on the OpenStreetMap PostGIS wiki page, they just install the postgresql-9.1-postgis package.  Yet this Debian page seems to imply that you should install both packages as postgis contains common PostGIS files while postgresql-9.1-postgis contains files that specifically support PostgreSQL 9.1

Thanks.


texts: Defining which layer (overlay) is displayed on Leaflet map?

texts:
I have 7 layers overlay on a base map, i need to define which layer is selected (active,displayed,shown) on the map.

var overlays = {
        "AL HOUDA Espaces_Verts": alhouda6_layer,
        "AL_HOUDA_imm": alhouda7_layer,
        "AL_HOUDA_villa": alhouda8_layer,
        "AL HOUDA_ECO": alhouda5_layer,
        "AL_HOUDA_EQUIP": alhouda4_layer,
        "AL_HOUDA_REC": alhouda3_layer,
        "AL_HOUDA_PARK": alhouda2_layer,
        "AL_HOUDA_200lOG": alhouda1_layer,
    };

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

For example when i select "AL_HOUDA_villa" and it's shown on the map i want to define that it's selected so i can show data just for this layer


texts: Determining if coordinates from feature are in EPSG:4326


texts: How to drag a point and then use savestrategy in OpenLayers

texts:
Im trying to use the wfs-t in openayers using the savestrategy.
I have created a savestrategy and it works when I for example uses the DrawFeature Control to add a point. A new point is added to the database.
However, when I use the DragFeature control to move a point (in the same vector layer as I have added using DrawFeature) it doesnt work.

The response from Geoserver is

<wfs:InsertResult>
   <ogc:FeatureId fid="none"/>
</wfs:InsertResult>
<wfs:TransactionResult>
  <wfs:Status>
    <wfs:SUCCESS/>
  </wfs:Status>
</wfs:TransactionResult>

And to me that looks like for some reasons no information sent regarding the fid that should have been moved.

What could I havce done wrong?
BR
Mike


texts: can the Identify Features button in QGIS list data in natural order

texts: Is there anyway to get the Identify Features button in QGIS to list my attribute data in their natural order, NOT in alphabetical order?


texts: How to extract a single road from OSM GIS Database?

texts:
Is there a way to extract a single row out of the OSM-GIS-Database (by a Query)?
It would be so much easier if I could find an EER Model of the Database.

I was searching a long term and found very less information about that. So I think it would be the best to ask the pros before thinking it's impossible.


texts:
When I am importing RASTERIO then i am getting an error as:

Traceback (most recent call last):
  File "<pyshell#0>", line 1, in <module>
    import rasterio
  File "C:\Python27\lib\site-packages\rasterio\__init__.py", line 17, in <module>
    from rasterio._base import eval_window, window_shape, window_index
ImportError: DLL load failed: The specified module could not be found.

On importing FIONA i am getting error as:

Traceback (most recent call last):
  File "<pyshell#1>", line 1, in <module>
    import fiona
  File "C:\Python27\lib\site-packages\fiona\__init__.py", line 78, in <module>
    from fiona.collection import Collection, BytesCollection, vsi_path
  File "C:\Python27\lib\site-packages\fiona\collection.py", line 7, in <module>
    from fiona.ogrext import Iterator, ItemsIterator, KeysIterator
ImportError: DLL load failed: The specified module could not be found.

While both rasterio and fiona installed successfully by pip install. please help some one i am stucked here.i am using python 2.7 and gdal1.10 and numpy1.10.2 on 32 bit


texts: Creating separate dxf file for each ArcSDE featureclass with FME 2017?

texts:
I am having problems creating separate dxf files for my sde Featureclasses with fme2017. I am not a CAD user but have to export over 200 featureclasses out of my 10.3.1 SDE Database. I have created the sde reader which splits the database into the featureclasses and would now like to create a separate dxf for each individual featureclass. Because I have over 200 I want to do this automatically. I have already tried using a featureclassfilter but the output is always a single dxf. If I have tried using a generic writer or  fan-out with featuretype as parameter it still routes all to the same single dxf.

How can I create one dxf for each Featureclass?


texts: Change feature layer symbol of single polygon on mouse hover

texts:
My goal is to display a complete feature layer of ~400 features with a default symbology when the map loads and then replace the symbology of a single polygon in that feature layer on a mouse hover event. Eventually, I'll want to be able to select that polygon and permanently alter the symbology while also changing the symbology of the remaining polygons in the same feature layer... this is basically a highlight and whiteout.

Here how I started:

map = new Map("mapDiv", {
  basemap: "topo",
  center: [-50, 50]
});

map.on("load", function(){
    map.graphics.enableMouseEvents();
    map.graphics.on("mouse-out", function(evt){
        map.graphics.clear();
    });
});

var layer = new FeatureLayer("http://myserver/arcgis/rest/services/my_map/MapServer/1", {
      mode: FeatureLayer.MODE_SNAPSHOT
});
map.addLayer(layer);

var hover_symbol = new SimpleFillSymbol(
  SimpleFillSymbol.STYLE_SOLID,
  new SimpleLineSymbol(
    SimpleLineSymbol.STYLE_SOLID,
    new Color([38,115,0]),
    2
  ),
  new Color([38,115,0,0.75])
);

layer.on("mouse-over", function(evt){
    var hover_graphic = new Graphic(evt.graphic.geometry, hover_symbol);
    map.graphics.add(hover_graphic);
});

This works to get a highlight over the feature that the mouse is hovering on. I now realize, however, that the new graphic is merely covering the old feature layer graphic instead of replacing it.

I see a Unique Value Renderer that allows you to style a feature layer by an attribute, but it seems like I would need to query the layer for that to apply here.

Any ideas on how to accomplish this on-the-fly? Can I somehow loop through the graphics in the feature layer and single out the polygon over which the mouse is hovering?


texts: For remote debugging of an ArcGIS Server extension, how can I identify which ArcSOC process is handling my target service so I can attach to it?

texts:
I'm remote debugging a ArcGIS Server 10.6 extension with Visual Studio.  I've got it set up so it all works but I'm having trouble identifying which ArcSOC process I should attach to.

I found this discussion:
How can I determine which service is associated with a particular ArcSOC.exe process and vice versa?

But that is for a much older version and doesn't apply as far as I can see.


texts: Brinkhoff Generator Esri street shapefile + nodes and edges files

texts:
This is the original Brinkhoff generator.: http://iapg.jade-hs.de/personen/brinkhoff/generator/

& this is a modified version for higher workloads: http://moto.sourceforge.net/

The modified version does two things as their website says:
1) It takes an ESRI shapefile and produces node- and edge files in the format of the Brinkhoff generator.
2) It takes node- and edge files and generates position and query traces.

I tried esri shape files with the node/edge generator and it does not seem to read it, and the node/edge files provided in the original Brinkhoff website does not work with this version as it includes nodes names.

I get this message: Exception in thread "main" java.lang.RuntimeException: No names allowed in nodes!

Does anyone know what esri shape files I can use (and if possible any free downloads?)

How can I remove nodes/edges names from the Brinkhoff  node/edge files?


texts: merging geometry (Feature Collections from Earth Engine) and Fusion table data for Supervised classification

texts:
I have data for a crop in fusion table( I had in excel format hence I converted to fusion table to use in Earth Engine). I also have feature collections (land cover classes like urban, plantation etc) obtained from earth engine platform digitization.
Now I need to merge the crop data and other land cover data for classification for giving training set  While I try to do classification, it is showing error :
Dictionary (Error)
Property 'lulc' of feature '2_2_0' is missing.
Please Give suggestions


texts: Creating markers perpendicular to different lines in QGIS

texts:
I want to create markers on different lines and this markers should have always 90° on the line.

I used the simple marker, with a 90° angle, but this doesn't work for every line.

How can I get a 90°angle for every line segment?


texts: Problem when compiling QGIS (with customizations) for Windows

texts:
I was able to customize a few features of QGIS in Ubuntu. Now I am trying to compile the same code on Windows, by following the installation guide given for Windows. I have downloaded all of the dependencies using OsGeo4W and struck at SIP. I tried setting environment variables to point to sip site-pacakges too ( C:\OSGeo4W\apps\Python27\lib\site-packages )
Any guidance would be of lot of help. I am using visual studio compiler libraries.

Here is the cmake trace

Quantum GIS version: 1.8.0 Lisboa (10800)
Could not find GRASS
Touch support disabled
Found Proj: C:/OSGeo4W/bin/proj.dll
Found GEOS: C:/OSGeo4W/bin/geos_c.dll
Found GDAL: C:/OSGeo4W/bin/gdal18.dll (1.9.2)
Found Spatialindex: C:/OSGeo4W/bin/spatialindex1.dll
Found Qwt: C:/OSGeo4W/bin/qwt5.dll (5.2.1)
Found Sqlite3: C:/OSGeo4W/bin/sqlite3.dll
Found PostgreSQL: C:/OSGeo4W/lib/libpq.lib
Found SpatiaLite: C:/OSGeo4W/bin/spatialite.dll
Pedantic compiler settings enabled
Found Python executable: C:/Python27/python.exe
Found Python version: 2.7
Found Python library: C:/Python27/libs/python27.lib
CMake Error at cmake/FindSIP.cmake:54 (MESSAGE):
  Could not find SIP
Call Stack (most recent call first):
  CMakeLists.txt:476 (FIND_PACKAGE)

Configuring incomplete, errors occurred!

Thank you.
Srikanth


texts: Is there a 1.5:1 ratio coordinate system?

texts:
I'm a student fairly new to using GIS. I am doing a college project on using ArcGIS to map virtual worlds (such as the continents in World of Warcraft). The game programmers used a custom coordinate system. "Map coordinates are a percentage of the width and height of these images" to quote this site: WoWWiki: Coordinate-safe maps

The maps are all rectangular images with the same proportion ( 1.5:1 ) The in-game distances and geography seem to match the visual appearance of the maps. It's a "flat earth" world rather than a sphere. So there is no visual distortion as there would be on an Earth projection. I want to use a coordinate system that matches the x-y coordinates used by the game and displays correctly on the map image used as a layer - without turning the image from an oblong into a square.

Is there a way to use or make a coordinate system in ArcGIS where the X coordinates are "stretched" relative to the Y ?

Tried to give as much detail as possible. Hope it makes sense.


texts: Cropping raster using multiple shapefiles

texts:
I have a raster which I want to crop multiple times using different shapefiles. My code uses mask method which returns numpy array and its metadata:

out_img, out_transform = mask(raster=data, shapes=coords, crop=True)

I run different kinds of thresholding on out_img and I want to crop it again using another shapefile but I'm not sure how to do that since out_img is simply a numpy array.

Do I need to save out_img and out_transform as TIFF, load it and use mask method again? Is there any more straightforward method of achieving it (i.e. without saving the intermediate steps)?


texts: Setting shapefile points size to actual acerage using QGIS?

texts:
I have a shapefile that has lat/long points data of specific facilities in the US. This data is joined with census blocks group population data because eventually I want to calculate how many people live within a 1 mile and 3 mile radius of each facility.

The issue I am having is I do not want just a uniform point for each facility. I want each point to be represented by how many acres the facility is (I have acreage data in the table).

Now in QGIS, I can arrange the size of each point based on acreage, but its not geographically accurate. It just shows that small acre plots are 1 pixel in width, and big acre lots are >1 pixel in width. (I did this by using the "Assistant | Symbol Size" tool in the "Symbols" tab of "Layer Properties"

How can I get each point to physically be the acre size, and then create the 1 and 3 mile buffers around that? So if there is a plot of land that is 25 acres, I want a circle that is exactly 25 acres in size. THEN place a 1 mile and 3 mile boundary on that.

Some photos to help explain what is going on.

Here, the points are all uniform (with the 1 mile boundary in purple, 3 mile boundary in blue), but I would like them to be based on the actual acreage. Like if the plot of land was 25 acres, I want a circle that is exactly 25 acres in size.

I can adjust the points like this, but its not physically and geographically real to the size of the acre. Now I understand that maybe the boundary of these plots of land are not uniformly circular, but that is an assumption I am making to make the analysis possible.

My methodology. need help with step 2

I am using QGIS.

EDIT


texts: Shapefile import used mapbuilder to oracle and publish on Geoserver but can't read '_'(underline)

texts:
environment:windows10 and oracle11g
tools:oracle mapbuilder11 and geoserver2.8.2
step:
1.use mapbuilder to shp file into oracle .  like this http://www.duoluodeyu.com/37.html
2.create a link with oracle like this http://docs.geoserver.org/2.5.x/en/user/data/database/oracle.html
3.geoserver add this to layer and publish
when i used mapbuilder to shp file into oracle. if the geomatry table name is not have '_'(underline)

everything is OK

But if the geomatry table name add '_' like XINZHUANG_2004
geoserver publish result is error ,like this

i don't know what is going on and I don't ready think it can't add '_'(underline).
I am so sorry for my poor english.


texts:
I am developing a mining sampling database.  There will be many project specific, user created projections.  Typically they will be Oblique Mercator defining a local coordinate projection system.

Does anyone have suggestions on managing a large number of custom SRIDS in PostGIS?

The easiest is to just add it to the public.spatial_ref_sys table but how will that complicate updates to PostGIS/PROJ versions?  Also what would be the best large integers to start my custom SRIDS to avoid future conflicts?

Or would it be best to create a public.spatial_ref_sys_custom table?  I am not sure that would be ideal since it would create numerous issues with on the fly st_transform() queries.

Has anyone else ran into this issue before?


texts: Is it necessary to calibrate LandSat-8 (L1T) for the calculation of TOA radiance prior to atmospheric correction, or is it already calibrated?


texts:
When I pipe the output of raster2pgsql to psql -d ais -U postgres -h localhost -p 5432 it asks for a password.

How can I supply it the password on the command line?


texts:
I have created a hexagonal grid shapefile using MMQGIS and want to use this for hexagonal binning.

I have another point shapefile which contains patient records for three different types of hospital activity (let's call them A,B and C) for two different years.

What would be the best way of quickly generating point counts (via Points in Polygon) for each subset of the activity shapefile so that I ended up with 9 additional columns in my hexagonal grid shapefile?

Year 1, activity A
Year 2, activity A
Difference between years 1 and 2 for activity A
Year 1, activity B
Year 2, activity B
Difference between years 1 and 2 for activity B
Year 1, activity C
Year 2, activity C
Difference between years 1 and 2 for activity C

Is it possible to script this using pyqgis? I've started reading up on pyqgis but I'm a long way off being able to tell if this is possible.

If it is possible, are there any tutorials out there that would cover this sort of processing or at least point me in the right direction to start figuring it out?

I really want to learn scripting in QGIS and this feels like something which would lend itself to getting started.


texts: How to access the options of the Label Properties dialog in QGIS2.14.2?

texts:
The Change Label tool has a Label Properties dialog but the options are greyed out.

How to access them?


texts: Making list of MXDs from multi-value parameter in ArcPy?

texts:
I would like to selectively choose my MXD documents through the multivalue or filter parameter in the toolbox.

From the selected MXD's the user should be able to replace text and re-save the document.

I can get it to work on a series of unselected MXD documents located in a specific folder, but I need to be able to choose those MXD's.

Essentially the question is how do I integrate the multi-value parameter into the code?

import arcpy
import os

inMXD = arcpy.GetParameterAsText(0) #toolbox parameter workspace
#prepare for multivalue toolbox parameter
#ListInMXD = inMXD.split(";")
#set old and new text variables
OldText = arcpy.GetParameterAsText(1)
NewText = arcpy.GetParameterAsText(2)

#THis is where it should loop through the ListInMXD

for file in ListInMXD:
    filePath = os.path.join(ListInMXD, file)
    mxd = arcpy.mapping.MapDocument(filePath)
    for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
        if elm.text == OldText:
            elm.text = NewText
            mxd.save()

With PolyGeo's answer I was able successfully to re-write the code as follows:

import arcpy
import os
mxdString = arcpy.GetParameterAsText(0) #toolbox parameter workspace
#prepare for multivalue toolbox parameter
mxdList = mxdString.split(";")
arcpy.AddMessage(mxdList)
#set old and new text variables
OldText = arcpy.GetParameterAsText(1)
NewText = arcpy.GetParameterAsText(2)
#THis is where it should loop through the ListInMXD
for mxdItem in mxdList:
    mxd = arcpy.mapping.MapDocument(mxdItem.strip("'")) #this little snippet came from https://geonet.esri.com/thread/67229
    for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
        if elm.text == OldText:
            elm.text = NewText
    mxd.save()
del mxd


texts: ArcSDE - use native geometry types for later migration away from SDE?

texts:
Background: We have a non-spatial reporting web application with a SQL Server backend.  We would like to add a spatial component to this application: serve up some basemaps, serve up points and polygons for view and editing, build reports for selected polygons.

To that end, the plan is to use ArcSDE with SQL Server, and use ArcServer to serve up base layers and vector geometries.  (something like WMS/WFS-T)

Problem: We may, at some point, decide we want to move away from expensive proprietary software, and use a FOSS alternative (e.g. GeoServer, or TinyOWS if/when it supports SQL Server).  If possible, I'd like to avoid building the application irrevocably around ESRI technology.

My understanding of SDE is that it can either use ESRI specific spatial formats, or spatial formats native to the database system it's installed on.  (I presume there is some performance penalty for using the native formats.)

My Question(s): If we were to set up SDE on this SQL server database using native SQL Server spatial formats, would we be able to set up say, GeoServer, to run against the tables that SDE was managing?  What would it take to rip SDE out and use something else?  (Also, a meta-question: is this a reasonable question to be asking?  Is there something important I'm missing?)


texts: How to convert shapely geometry linestring into shapefile?

texts:
I have extracted linestring from Landsat images. now, I am literally stuck what to do? I want to convert it into shapefile. It is actually coastline geometry.
   **In[1]:** type(water_lines_smooth)
    **Out[1]:** list
    **In[2]:** print(water_lines_smooth)
    **Out[2]:**
    [<shapely.geometry.linestring.LineString at 0x2d11cd01c08>,
     <shapely.geometry.linestring.LineString at 0x2d11019f048>,
     <shapely.geometry.linestring.LineString at 0x2d11cd09d08>,
     <shapely.geometry.linestring.LineString at 0x2d11cceac08>,
     <shapely.geometry.linestring.LineString at 0x2d11ca7d0c8>,
     <shapely.geometry.linestring.LineString at 0x2d11ca7d108>,
     <shapely.geometry.linestring.LineString at 0x2d142747f88>,
     <shapely.geometry.linestring.LineString at 0x2d142747ec8>,
     <shapely.geometry.linestring.LineString at 0x2d142747088>,
     <shapely.geometry.linestring.LineString at 0x2d11ca7d148>,
     <shapely.geometry.linestring.LineString at 0x2d142747e48>,
     <shapely.geometry.linestring.LineString at 0x2d142747d08>,
     <shapely.geometry.linestring.LineString at 0x2d142747d48>,
     <shapely.geometry.linestring.LineString at 0x2d142747cc8>,
     <shapely.geometry.linestring.LineString at 0x2d142747408>,
     <shapely.geometry.linestring.LineString at 0x2d142747e08>,...

    **In[3]:** [{'geometry': x.__geo_interface__} for x in water_lines_smooth]
    **Out[3]:**
    {'geometry': {'type': 'LineString',
       'coordinates': ((377205.43543222273, 2464301.85113062),
        (377235.43946816603, 2464301.85113062),
        (377235.43946816603, 2464241.855100152),
        (377205.43543222273, 2464241.855100152),
    (375315.1811677954, 2464391.8451763224),
        (375285.1771318521, 2464331.8491458544),
        (375225.1690599655, 2464331.8491458544),
        (375165.1609880789, 2464361.8471610886),....
I have tried to convert it into geodataframe but there is geometry column issue.

     **In[4]:** df=gpd.GeoDataFrame(water_lines_smooth)
     I**n[5]:** df.to_file("shapef.shp")
    AttributeError: No geometry data set yet (expected in column 'geometry'.) –

is there any alternative way to get shapefile?


texts:
I have a .tif raster file which contains integer values at 1, 2, 3, 4 and 5. Importing this file into ArcMap also shows a number of cells with Nodata value where the above classes are not present.
I was able to convert a tile from this .tif file into a numpy array using the Arcpy documentation. I basically ran the following code
point = arcpy.Point(XCoord, YCoord)
test = arcpy.RasterToNumpyArray(RasterFileName, point, 10,10)

Where XCoord and YCoord are specifically selected to be in an area where all the returned raster values should be Nodata (according to ArcMap). Instead, I get the following output when I print out test:
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

Why are these not the values corresponding to ArcMap?. For another set of coordinates in the region where there is data I get the following:
array([[2, 3, 2, 1, 2, 2, 2, 1, 1, 2],
       [3, 3, 2, 1, 2, 2, 2, 1, 2, 2],
       [3, 2, 2, 2, 2, 2, 2, 2, 3, 4],
       [3, 2, 2, 3, 2, 3, 3, 3, 3, 4],
       [2, 2, 3, 3, 2, 2, 3, 2, 2, 2],
       [1, 3, 3, 3, 2, 2, 2, 1, 1, 2],
       [1, 2, 2, 2, 2, 2, 2, 2, 1, 2],
       [1, 2, 2, 2, 2, 2, 3, 3, 2, 2],
       [1, 1, 1, 1, 1, 2, 4, 4, 2, 2],
       [1, 1, 1, 1, 1, 2, 3, 3, 3, 3]])

This corresponds to the ArcMap information in that tile. Therefore I am sure that the raster is converting correctly.
How does ArcMap differentiate between NoDATA and 0 valued raster cells when in numpy NoDATA values are not allowed? Is the returned numpy array the actual data or is it simply a higher level representation of something else?


texts: Is there any specific index for identifiying settlements in remote areas through Landsat or Sentinel images?

texts: I am trying to identify (especially) small settlements in Nigeria that have no electrification. So is there any index that could help me detect households from Landsat or Sentinel images?


texts: Redirect remote shell command output to ArcMap geoprocessing tool display?

texts: I have a Python script that runs a shell command on a remote server (using the Fabric module).  When I run the tool in the debugger, I see all the output scrolling, as the command can take hours to run.  I have created an ArcMap tool to run from the ArcToolbox, but I cannot see the progression of the shell command's output.  The geoprocessing window just stays open until finished.  Is there anyway to send the shell output to the geoprocessing window in real time?


texts: Unable to cut and paste between 2D & 3D layers in QGIS

texts:
I have two point layers, one with a subset of the fields of the other.  I want to cut and paste points from the layer with fewer fields into the other but QGIS will not allow me.

Source layer is created in memory by a python script and has a 2D geometry.
Destination layer is a postgres table with 3D geometry.

I get a banner message across the top of the window saying No points copies and the message vanishes before I can study it in detail.  There are some symbols that might be buttons but I can't tell.

The issue is that I can't figure out what is wrong as there are no error messages that I can find in the logs.

How can I diagnose the issue?

QGIS 2.18.7 on a Mac.


texts: Really slow extraction from raster even after using crop?

texts:
I have a large raster file (245295396) cells and stacks of rasters having 4 layers each which lie in the extent of this large raster. To start with I am trying to get value from one stack (3 channels) and for the same zone from the large raster. Every things works fine, just the extraction from large raster takes 5 mins. So, if I repeat this process for 4000 more times it will take 13 days.

cld<- raster("cdl_30m_r_il_2014_albers.tif") #this is the large raster
r<- stack(paste(path,"/data_robin/", fl,sep="")) #1 stack,I have 4000 similar
mat<-as.data.frame(getValues(r)) # getting values from the stack
xy<-xyFromCell(r,c(1:ncell(r)),spatial = TRUE)
clip1 <- crop(cld, extent(r)) # Tried to crop it to a smaller size
cells<-cellFromXY(clip1,xy)
mat$landuse<- NA
# mat$landuse<-cld[cells]
mat$landuse<- extract(clip1,cells) #this line takes 5 mins based on profiling

> cld
class       : RasterLayer
dimensions  : 20862, 11758, 245295396  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : 378585, 731325, 1569045, 2194905  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
data source : /Users/kaswani/R/Image/cdl_30m_r_il_2014_albers.tif
names       : cdl_30m_r_il_2014_albers
values      : 0, 255  (min, max)

> r
class       : RasterStack
dimensions  : 9230, 7502, 69243460, 4  (nrow, ncol, ncell, nlayers)
resolution  : 0.7995722, 0.7995722  (x, y)
extent      : 589084.4, 595082.8, 1564504, 1571884  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
names       : m_3608906_ne_16.1, m_3608906_ne_16.2, m_3608906_ne_16.3, m_3608906_ne_16.4
min values  :                 0,                 0,                 0,                 0
max values  :               255,               255,               255,               255

My data is in .tiff format and I am new to geospatial coding.

I have also tried the approach at Increasing speed of crop, mask, & extract raster by many polygons in R? but during the masking part it gives an error Error in compareRaster(x, mask) : different extent.


texts:
I have an Mapserver 7.0.4 for WMS GetMap requests currently over GET.
Is it possible to send POST GetMap requests to the mapserver?
I already searched some other topics and the documentation, it seems that WFS supports POST but i have no information about WMS.

I already tried to send some POST requests but all i got back was the template i specified and not the actual image i expected.
(Capabilities says POST is enabled)

But it seems that the Mapserver is indeed processing the POST Body:
In the POST Body i used a modified XML file i found online, if the xml file contains mistakes it gives out an error. So something is definitely going on there.

GET: example URL (which works fine)

http://exampleurl.com/cgi-bin/mapserv?MAP=/etc/mapserver/example.map&LAYERS=Example&SRS=EPSG:4326&VISIBILITY=true&TRANSPARENT=TRUE&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&FORMAT=image/png&BBOX=7,49,8,50&WIDTH=256&HEIGHT=256

POST: XML (Which does not work)

<?xml version="1.0" encoding="utf-8"?>
<ogc:GetMap gml="http://www.opengis.net/gml"
 xmlns:ogc="http://www.opengis.net/ows"
 xmlns:gml="http://www.opengis.net/gml" version="1.3.0" service="WMS">
<StyledLayerDescriptor version="1.0.0">
<NamedLayer>
  <Name>Example</Name>
</NamedLayer>
</StyledLayerDescriptor>
<SRS>EPSG:4326</SRS>
<BoundingBox srsName="http://www.opengis.net/gml/srs/epsg.xml#4326">
<gml:coord>
  <gml:X>7</gml:X>
  <gml:Y>49</gml:Y>
</gml:coord>
<gml:coord>
  <gml:X>8</gml:X>
  <gml:Y>50</gml:Y>
</gml:coord>
</BoundingBox>
<Output>
<Size>
  <Width>256</Width>
  <Height>256</Height>
</Size>
<Format>image/png</Format>
<Transparent>TRUE</Transparent>
</Output>
<Exceptions>application/vnd.ogc.se+xml</Exceptions>
</ogc:GetMap>

If WMS works with POST, has anyone more information about the structure of these XML files and how to build them properly?
The Mapserver documention contains only:

  The ability
  for MapServer to be able to receive Post requests with XML-encoded information sent in the body of the request has
  been added. Also, the ability to generate XML-encoded Post requests for WFS layers has been added.


texts: In the Symbology Properties of an ArcGIS layer, why is the count sometimes not shown?

texts:
Sometimes the count column does populate for me, but most times it does not seem to. I haven't been able to figure out a rhyme or reason to it. There is no definition query or anything out of the ordinary. Is this a bug or is there something you can do to make sure that field is populated?

I am running ArcView 10 on Windows 7 x64.


texts: WMS feed displays in ArcCatalog 9.3 but not in ArcMap 9.3

texts:
I am facing a very very weird behaviour. We used to work with a WMS feed in ArcGIS 9.3 and, since a few weeks, we are no longer able to display it.

The providers of the feed informed us that they have "made some changes" in the feed, which might cause the problem. They have now reverted some of the changes and the behavior is as follows:

ArcCatalog 9.3: the WMS sublayers can be displayed
ArcMap 9.3: the WMS sublayers appear in the TOC but the data view remains "white and empty" (no error message is displayed)
ArcCatalog 10.1: the WMS sublayers can be displayed
ArcMap 10.1: the WMS sublayers can be displayed

As anyone been facing something like this? Do you have any idea how I can monitor it or what can cause this incompatibility in 9.3? I believe that the server should send back the WMS tiles along with an XML file. Anyone knows where this XML file would be stored by ArcGIS?

The WMS feed is served with MapServer 6.3

Thanks a lot in advance!


texts: How can we sync attribute tables with excel file automatically?

texts: I always export to excel every time I update a shapefile, is there a way of updating excel automatically.


texts: Mean lower low water (MLLW) for tidal correction

texts: I am estimating water depth (bathymetry) using satellite data over a US water. I have obtained tidal values (datum is MLLW) as at the time of satellite overpass to correct my satellite-derived bathymetry for tide. I was wondering if the value obtained (e.g. 0.68m) will be subtracted or added to my SDB estimates - to account for tide.


texts:
I have a long list of GeoHash codes for a certain geographic region. I'm using 6-character precision so for example if the area under study was Earl's Court in London, I would have a few six letter codes that for that region.

What I want to do is to be able to generate heat maps using this information. Specifically I want to colour a set of codes based on a certain condition.

I have searched but the only thing I have seen is the GeoHash query page (http://geohash.org/c216ne), but that just returns a marker to the a specific position, not the area.

Are there any tools for generating a heat map with this geoHash Information?


texts: Create geometry from point data - QGIS 3.4.1

texts:
Is there an efficient way to produce geometry from point data using QGIS?

Each point represents a feature that may require one or more geometric primitives such as circles, curves or polylines (open or closed) to represent it.

Please note that while CAD is the best tool for creating geometry I am seeking alternative methods of achieving this using QGIS without using CAD. However the final QGIS geometry needs to be exportable to CAD in a usable format.

The goal is to produce database linked maps within QGIS but be capable of supplying useful CAD geometry to clients if they request it.

Options explored so far:
 Vector > Geoprocessing Tools > Buffer
This option is a good start that can approximate circles, however the geometry is segmented and when exported to CAD (using DXF) it is interpreted as hatching. Hatching requires additional processing within CAD to produce usable geometry and being restricted to segmented circles also creates problems. This approach is less efficient than programmatically creating geometry within CAD from raw point data.

Python scripting
I'm barely scratching the surface here as there is a lot to learn. If this is the best option then any advice on how to go about it is very helpful.

Edit: This question has been rewritten to clarify the goals and outline what has been explored.


texts:
Yesterday ArcMap crashed and corrupted the toolbar layout of a map I used as model.

No TOC, the Menu bar is on the left, ... Nothing dramatic, so I replaced them to something that is more of my liking:

I saved, closed my map (to save the display in my Normal.mxt) and reopened the map. I was back at its corrupted "look"!

I deleted my Normal.mxt template. Open ArcMap, and it was ok. I added the toolbars I use the most, saved it and made a copy of my Normal.mxt.

If I open an new map, the display is fine. If I try to open any existing map (even one created a couple months ago), the display is "re-corrupted" but I know I was fine in that document.

I used a script to create 45 maps from the model mentioned earlier, so I expect to them to have the corrupted layout (didn't notice this before creating them). But I don't want to have to replace all my toolbars everytime I open a map I created before.

Never had that behavior before.

I look like it's the same question as "ArcMap changing toolbars", but I didn't find an answer there for me.


texts: Converted GPX file has null values in open attributes

texts:
I'm an absolute newbie but I think I'm doing this right.

I've got a bunch of .fit files from my garmin. I can use GPSBabel to manually convert these to GPX XML and the map looks fine after the conversion. If I open the file it looks fine.

I then use Vector | GPS | GPS Tools to import the gpx file. The layer is created but there's nothing to see on the map. If I go to the Attribute Table there is only one line and it's full of nulls.

Is there are better output file type in GPS Babel that qgis prefers?

<?xml version="1.0" encoding="UTF-8"?>
<gpx version="1.1" creator="GPSBabel - http://www.gpsbabel.org" xmlns="http://www.topografix.com/GPX/1/1" xmlns:gpxx="http://www.garmin.com/xmlschemas/GpxExtensions/v3" xmlns:gpxtpx="http://www.garmin.com/xmlschemas/TrackPointExtension/v1">
  <metadata>
    <time>2014-04-25T05:42:52.371Z</time>
    <bounds minlat="-36.755794164" minlon="174.555711483" maxlat="-36.722196898" maxlon="174.580664325"/>
  </metadata>
  <trk>
    <trkseg>
      <trkpt lat="-36.755794164" lon="174.563342619">
        <ele>42.200000</ele>
        <time>2014-04-24T22:28:41Z</time>
        <extensions>
          <gpxtpx:TrackPointExtension>
            <gpxtpx:atemp>27.000000</gpxtpx:atemp>
          </gpxtpx:TrackPointExtension>
        </extensions>
      </trkpt>
      <trkpt lat="-36.755728199" lon="174.563325436">
        <ele>42.200000</ele>
        <time>2014-04-24T22:28:45Z</time>
        <extensions>
          <gpxtpx:TrackPointExtension>
            <gpxtpx:atemp>27.000000</gpxtpx:atemp>
          </gpxtpx:TrackPointExtension>
        </extensions>
      </trkpt>


texts: How do I run an extended overpass turbo query from a script?

texts:
I have an xml query which starts out like this:

<osm-script output="xml" timeout="250">
  <id-query {{nominatimArea:Ridgefield Connecticut}} into="area"/>
  <!-- gather results -->
  <union>

This works fine here: http://overpass-turbo.eu/s/5J7

But when I try to query using wget, it returns HTTP request sent, awaiting response... 400 Bad Request. If I resolve the query into it normalized form, then it works fine:

<osm-script output="xml" timeout="250">
  <id-query type="area" ref="2705805056" into="area"/>
  <!-- gather results -->
  <union>

How can I pass the former query through a script, or how can I lookup the area in advance?


texts: Converting GPS track to turn by turn directions

texts:
Working for a client on a tracking package for delivery vehicles, I've been given an interesting requirement.  Over the past several months, we've been collecting telemetry from the vehicles, including lat/lon, speed, heading, and so forth.  The client has asked if we can take an arbitrary track (defined here as an ordered set of lat/lon coordinates) and convert it to a series of turn-by-turn instructions that would result in the same route.

I've been fiddling with OSRM and Graphhopper and both are able to pretty easily turn the points into a polyline and return a set of instructions, but the instructions need a serious amount of reduction in order to be useful to a driver -- for example, traveling down a street for half a mile may result in a dozen instructions that just say "continue on _______ street for 200ft" or whatever the interval between the data points happens to be.

Before I just "brute force" the answer by removing all intermediate points where the street name is the same, I am trying to find an existing SDK or library that can do this for me.  My concern is that blindly removing points may result in generated directions that take a different path, and in this situation even a slightly different path is unacceptable.  At present I'm not even sure what I would call a feature like this, and all I've been able to find when searching are websites that can do this for you, manually, for cyclists and the like.

Edit: Example lat/lon pairs

43.059595, -70.753773
43.060723, -70.753812
43.061808, -70.754822
43.062661, -70.755752
43.063645, -70.756604
43.064411, -70.756960
43.065192, -70.757475
43.065901, -70.757851
43.065858, -70.758326
43.065597, -70.758999
43.065250, -70.759751
43.064888, -70.760860
43.064700, -70.761494
43.065337, -70.761810
43.065988, -70.762088
43.066581, -70.762345
43.067072, -70.762563
43.067376, -70.762701
43.067463, -70.762642
43.067521, -70.762206
43.067651, -70.761632
43.067752, -70.761276
43.067868, -70.760781
43.068027, -70.760444
43.068186, -70.760009
43.068403, -70.759593

Given to /route results in this polyline:

mayeGpaznLaFVwEjEiDzDgE`DqCxAaDxAmCfABfBv@`C|@zCfAzEl@xB}B`AaCv@wBn@aBh@{@XUGOsAYqBOiAWaB[eAc@qAg@uA

From this json result: https://pastebin.com/qLHrxVPu

Which results in many repeated instructions to continue forward on the same road, for example.  What I'm looking to create, from existing sets of lat/lon points like that, is a simple set of driving instructions saying "go straight on this road xxx meters, turn left here, turn right there..."

I know I can brute force this, but I'm hoping there's something built in somewhere.


texts: How does Geoserver compute bounding boxes from data?

texts: When creating a layer, Geoserver has the option of calculating bounding boxes from the data. I have only succeeded in doing this when I pass a Geometry type object in my SQL query. How can I define the SRS and bounding box when passing only, lon, lat, elevation, for example? Is this possible or do I have to use Geometry objects?


texts:
I am creating set of charts in arcmap using data driven pages
I can highlight the active page in the key map as shown in the picture but I want to create a dynamic leader to point to the current page. is that possible?


texts: How to create a webmap showing point density?

texts:
I am trying to develop hotspot map showing point density in our GIS based web application.

After searching for an option I have decided to use v.kernel command. I have a table containing some point data i.e. crime data. I need to display the hotspot of crime for a particular area.

In v.kernel vector file containing points and polygon are used for displaying hotspots/heatmaps. In QGIS there is an option for selecting some rows and exporting it to shapefile. Is there any option for running this command from a Python file or is there any better option for displaying hotspot/heatmaps?

I have tried OpenLayers heatmap, Python's heatmap and django-heat app but nothing seems to work for me. I am using OpenLayers TMS (Openlayers-2.12) for displaying maps generated through Mapnik and TileCache.


texts: Error while kriging within QGIS (R script)

texts:
I am trying to incorporate an R script within QGIS and have some problems with kriging interpolation. Here is my code :

G_STRUCT=gstat(NULL, id="OBS", formula=field~1, data=PointOrig, model=v_fit)
coordinates(PointOrig)=c('X','Y')
proj4string(PointOrig)=CRS(proj4string(PointGrid))
krige=predict(G_STRUCT, PointGrid)

NB : both PointOrig and PointGrid are point vector layers.

When I run the code, I get this error :

Error in bbox(dataLst[[1]]$data) : object not a >= 2-column array
Calls: predict ... predict.gstat -> getMaxDist -> apply -> bbox -> bbox

I tested the script under R Studio and it works well.


texts: How to print .dwg to .pdf with teigha true view

texts:
I don't have Autocad software and I can only find this free tool which can't read .dwg files.

Now, I need to convert .dwg to .pdf. But the teigha true view can't export .dwg

to .pdf.  But I find it offers a plot command. Unfortunately, I aslo don't have the

DWF6 ePlot.pc3 printer. So, does anyone know how to handle that or is that an impossible

task?  Do you know any other software can convert .dwg to .pdf?


texts: How to improve positioning of images taken with GPS and Gyro info?

texts:
I'm sort of 'recreating' google street view with MIT right now.

We have a car with a camera on top and we're driving at 10mph down a road. As we drive we're shooting pictures at 30fps and each photo is tagged with a timestamp. I can interpolate between images and GPS route so I can tell you the "estimated" position at which every image was taken.

However, we have some GPS errors, especially in downtown areas with multipath, etc. To take a first stab at this, we bought a better GPS unit from Trimble the Aardvark GPS + DR which is a decent GPS unit as well as a Gyro. We also tied in the wheel speed pulse to the unit so it can estimate the position of the car even if GPS drops out for a second or two. Over time though, drift occurs and it can tell us that we're 10-15 meters off the road until we get a decent GPS fix again.

We've done some basic metrics for the position, filtering out positions that have a faster than expected speed, given the previous points, and that can give me a decent guess at what parts of the route are bad.

That said, is there any way I can improve things further? I was going to do some map matching but I wonder if I should only do that in areas I decide are "bad" or if map matching is good for everything.

The real purpose of this is to be able to type in an address and pull up all images of that address. Right now we're close, but not close enough.

P.S.: Our current stack is Python, Google Maps and OSM mostly.


texts:
I am finding when I use label vector lines using the newer labeling tool that the labels always appear on top of the lines of that layer. And they appear on top of all other layers in my project. Is there a way to have these labels lower down (i.e. not on top of other layers).

Thanks,
David


texts: Using Search Cursor and Messaging Total Count of Items from Field in Layer that was Searched?

texts:
My program works but it is not displaying a message of the total sheet count found in attribute table that the search cursor searched through. For Example: It goes through the mxd chosen by user finds the layer with the name Grid in it and then is supposed to count the rows in the field Sheet_ID. End result is to print a message of the count found.  My code is below:

import arcpy, os, sys
from os import path as p
from arcpy import mapping as m

mxdList = arcpy.GetParameterAsText(0).split(';')

sheet_count = ()

for mapDoc in mxdList:
    arcpy.AddMessage(mapDoc)
    mxd = arcpy.mapping.MapDocument(mapDoc)
    for lyr in m.ListLayers(mxd, "*, Grid"):
        if lyr.description == "*, Grid":
        max_list = []
        rows = arcpy.SearchCursor(lyr)
        for row in rows:
            max_list.append(row.Sheet_ID)
        sheet_count = max(max_list)
        arcpy.AddMessage('Sheet count: %s'%sheet_count)
        print 'Sheet count: %s'%sheet_count
    else:
        arcpy.AddError('No Layers in %s match data source'%mapDoc)


texts: Wrong scale when importing dxf into QGIS?

texts:
Has anyone had the problem that the scale of the shp layer converted originally from a .dxf file turns out to be wrong?

For example, I draw a line of 10 meters long in AutoCAD and imported it in QGIS, however the length of this line turns out to be several thousand kilometres.


texts: Identify Number of Ground Control Points in DSM Using GDAL

texts:
I am trying to use GDAL (in Python) to extract metadata from a digital surface model (.tif).   I am using the "gdalinfo" command and applying different flags (e.g. -stats, -mm) to get additional information on the dataset, but I would like to know if it is possible to obtain the information on the number of ground control points used to help validate the DSM?  Per the GDAL documentation, this should be retrievable using the 'gdalinfo' command, but I am just not sure how to implement or if there are any additional flags I need to specify.
Is this actually feasible using GDAL, if control point metadata is actually embedded within the GeoTiff?


texts: Start/stop Python plugin during execution by just clicking on plugin button

texts:
I created a Python plugin in QGIS which is working fine. Now I want that whenever a user clicks on the plugin button (located on top of QGIS bar) it counts the number of clicks. If the number of clicks is odd my plugin should start running and if the number of clicks is even it should stop running (interrupting its execution). How can I do this and how can I map my code to this button so that it knows the number of clicks?

What I've done so far

    def initGui(self):
        icon_path = ':/plugins/my_tracker/icon.png'
        self.iconAction = self.add_action(
        icon_path,
        text=self.tr(u'My tracker plug'),
        callback=self.run,
        parent=self.iface.mainWindow())

    def run(self):
        if (self.iconAction.isChecked()):
            print ('Checked')
            self.show_markers()  // a method to do something

        else:
            print ('Unchecked')
            print("going to close socket", mySocket.close())

The problem is when I click on the plugin button it always runs the else part. It's not running the if part. Why?


texts:
Couldn't find any information about creating YSLD file (for styling my map, because it easiest way to add everything, what I want). I am working on the Geoserver site, and as "style" I can create only sld file.
How/where is it possible to create YSLD?


texts: Measuring QGIS performance to render 1-bit tif raster grids?

texts:
I notice that since weekly 2.17, it renders my 1-bit tif raster grids slower than 2.16, when I move viewport.
The new release 2.18 does the same as 2.17.
So I move back to 2.16.3.

Are there changes in the renderer that might cause this?

If yes, is there any how to configure to the old renderer?

Is this should be filed as a bug at http://hub.qgis.org/?

This bug has been fixed on 2016-11-4 0515 UTC by Martin Dobias.
The issue is closed now.

http://hub.qgis.org/issues/15796


texts: Incorrect color of raster data In GeoServer?

texts:
I publish a raster data with GeoServer, but a strange thing happened. If I publish the data with GeoTIFF, it is correct, like this:

But when I transform it as pyramid image and publish the data with ImagePyramid.The color is not correct. It become read like this:

When I zoom in the map, it became correct, like the GeoTIFF map.

There are another problem. Maybe the reason is the same as before. I have another raster data. I publish it as pyramid image, a black image in the map. If I zoom in the map, the image will appear.


texts: How to successfully do "nested" looping over a Fiona collection

texts:
I wrote some code that basically loops over all features in a shapefile with Fiona in Python, doing some things, calculating some things, etc.

I found that it would exit the for loop prematurely, depending on the data set I fed it. I found that it would work fine if I removed some code that (to me) has nothing to do with the loop we are in. This is my first attempt at GIS with Python, so I might be doing some newbie mistage with data structures I guess?

Below is a minimal example using Natural Earth data which you can grab at http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip

import fiona
from shapely.geometry import shape, Point

def find_containing_area(point, collection):
    for record in collection:
        geometry = shape(record['geometry'])
        if point.within(geometry):
            return record
    return None

with fiona.open('/tmp/ne_110m_admin_0_countries.shp', 'r') as source:
    for admin_record in source:
        properties = admin_record['properties']
        print(properties['name'])

        # in my original code, i created a point related to the admin_record
        point = Point(0,0)

        # if the following line exists, the loop will stop after the first
        # admin_record in this example. In my original code, it stopped
        # after just few records, not the first one though.
        containing_area = find_containing_area(point, source)

If I run it, it prints "Afghanistan" instead of the whole list of ~177 names. If I comment out the last line, it prints them all as intended.

My only guess is that someone an iterator is passed along with the collection? The only hint in the documentation that I found is "Fiona’s Collection is like a Python file, but is iterable for records rather than lines." which does help stupid-me not very much. http://toblerity.org/fiona/manual.html#data-model

edit: I have tried to copy the collection and pass the copy to the function instead, got the same bevahiour:

newsource = source[:]
containing_area = find_containing_area(point, newsource)


texts: GDAL c# API returned different WKT from gdalinfo

texts:
I've found completely different results in WKT strings while using gdalinfo utility and method from c# API GetProjectionRef(). In the first case the result was:

PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +
x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (4042799.739385183900000,6382762.252179093700000)

but in the second one

LOCAL_CS[\"WGS 84 / Pseudo-Mercator\",
    GEOGCS[\"WGS 84\",
        DATUM[\"unknown\",
            SPHEROID[\"unretrievable - using WGS84\",6378137,298.257223563]],
        PRIMEM[\"Greenwich\",0],
        UNIT[\"degree\",0.0174532925199433]],
            AUTHORITY[\"EPSG\",\"3857\"],
        UNIT[\"metre\",1]]

So, the questions are
1) Why results are different?
2)  How I can receive WKT  with projected coordinate system from c# code as in the first case?


texts: Shortest path (point to point) when filling out coordinates in boxes

texts:
I am using the processing toolbox interface to launch the "Shortest path (point to point)" command to connect two points using a line network layer (see attached photo).

If I fill out the starting and ending point coordinates of the Parameters Window by clicking directly in the screen, the algoritm works fine, but if I write out directly the coordinates of the two points in the Start Point and End Point boxes

I always get

There is no route from start point to end point.

Why I am getting this error if there is a route to connect both points?
I am using Windows 10, I have 16 GB of memory, Intel Core i7 vPro. QGIS version 3.12.0-Bucuresti


texts: Does PostGIS/pgrouting actually come with street info?

texts:
I've been looking at getting a new GIS setup going, involving QGIS/PostgreSQL/PostGIS/pgrouting, and so far it has been relatively easy. One thing I want to be able to do is to work out travel times, but before I go ahead and learn all the SQL I want to know, does PostGIS/pgrouting actually come with street info? It seems a bit optimistic that they would contain street directories for the entire world, so my question is where does one actually get that info? Are there open source databases available that contain that kind of information?

Also, I really only need Australian data, but being a data junky I won't say no to anything extra :D.


texts: H3-Python: How to generate hex bins over shapefile?

texts:
I am trying to generate hexbins over my shapefile to eventually cluster other geospatial events to them using H3. Here is my process, but I am wondering if there is a better way.
# packages
from h3 import h3
import pandas as pd
import geopandas as gpd
from shapely import geometry, ops
from shapely.geometry.polygon import Polygon
from shapely.geometry.multipolygon import MultiPolygon
from collections import Counter
from shapely.ops import cascaded_union

# transform multipolygon to polygon
def explode(indata):
    indf = gpd.GeoDataFrame.from_file(indata)
    outdf = gpd.GeoDataFrame(columns=indf.columns)
    for idx, row in indf.iterrows():
        if type(row.geometry) == Polygon:
            outdf = outdf.append(row,ignore_index=True)
        if type(row.geometry) == MultiPolygon:
            multdf = gpd.GeoDataFrame(columns=indf.columns)
            recs = len(row.geometry)
            multdf = multdf.append([row]*recs,ignore_index=True)
            for geom in range(recs):
                multdf.loc[geom,'geometry'] = row.geometry[geom]
            outdf = outdf.append(multdf,ignore_index=True)
    return outdf

# transform multipolygon to polygon
explode_vals = explode('C:\\Users\\Andrew\Desktop\\af_shape\\9c21fbef-c5ee-4eb6-88f3-ae4b938014bf202049-1-xml0al.1eimn.shp')

# set hex size
APERTURE_SIZE = 6

# Unify the boundries
union_poly = cascaded_union(explode_vals.geometry)

# Find the hexs within the shape boundary using PolyFill
hex_list=[]
for n,g in enumerate(union_poly):
    #print(n,'\r')
    temp  = geometry.mapping(g)
    temp['coordinates']=[[[j[1],j[0]] for j in i] for i in temp['coordinates']]
    hex_list.extend(h3.polyfill(temp,APERTURE_SIZE))

# create hex dataframe
hex_col = 'hex{}'.format(APERTURE_SIZE)
dfh = pd.DataFrame(hex_list,columns=[hex_col])
print('Sanity Check\nnumber of hexes:', len(hex_list))
print('number of duplicates:', len(hex_list) - len(dfh.drop_duplicates()))

# add lat & lng of center of hex
dfh['lat']=dfh[hex_col].apply(lambda x: h3.h3_to_geo(x)[0])
dfh['lng']=dfh[hex_col].apply(lambda x: h3.h3_to_geo(x)[1])

# create Point object based on hex latlng
dfh['geometry'] = dfh.apply(lambda x: geometry.Point(x.lng,x.lat),1)
dfh.crs = {"init": "epsg:4326"}

# plot hex latlng
dfh.plot(x='lng',y='lat',style='.',figsize=(48,48));

The second problem I have, is, if I want to filter my country level shape file to just a province, I get an error saying that my polygon is not iterable.
I am filtering like so, and following the same process.
explode_vals = explode_vals.loc[explode_vals['PROV_34_NA'] == 'Hilmand'].reset_index(drop=True)

I have also tried to use QGIS to export just the province of interest to geoJSON and a shapefile, but something about cutting it off from the rest of the country makes this whole process not work.


texts: Error using pip install for arcpy: "Could not find a version that satisfies the requirement arcpy...."

texts:
In my command line I typed

python -m pip install arcpy

to install arcpy. I got this error message in return:

"Could not find a version that satisfies the requirement arcpy (from versions: )
No matching distribution found for arcpy"

I'm using python 2.7. What to do?


texts: Show clustered features beside each other in OpenLayers

texts:
I have a map with clustered features, like in the picture. As you can see in the image some points have the same coordinates. What I want to do is when I'm at a certain zoom level and the points have identical coordinates I want to split them up so that the are beside each other so that the user can click on one of them.
Is this something that is possible to do with OpenLayers and clustered features?
const clusterSource = new olCluster({
            distance: 36,
            source: heatmapSource
});

const styleCache = {};
        this.clusterLayer = new olVectorLayer({
            visible: false,
            source: clusterSource,
            style: (feature) => {
                const size = feature.get('features').length;
                let style = styleCache[size];
                if (!style) {
                    style = new olStyle({
                        image: new olCircle({
                            radius: 8,
                            stroke: new olStroke({
                                color: '#000'
                            }),
                            fill: new olFill({
                                color: [153, 0, 0, Math.min(0.9 + (size / Math.max(positions.length)))]// '#3399CC'
                            })
                        }),
                        text: new olText({
                            text: size.toString(),
                            fill: new olFill({
                                color: '#fff'
                            })
                        })
                    });
                    styleCache[size] = style;
                }
                return style;
            }
        });


texts: Creating aggregate attribute from two polygon layers

texts:
I have two layers: one underlying land use information layer (with columns 'land_use_type' and 'area') and one for leased land.
To organize the leasing of the land it would be important that when drawing a polygon for an area I want to lease to have columns which tell me how much area of the drawn polygon is covered by land-use type (e.g. wood, farmland etc..).
I tried a bit around with the aggregate function and a virtual field, but to no avail:
aggregate(layer:='Land_Use',aggregate:='sum',expression:="area",filter:=contains(geometry(@parent), $geometry),order_by:="land_use_type")

How can I achieve that?


texts: openlayers. Get Data from feature when clicked (Modify events)

texts:
I have a folowing test case and this issue: https://github.com/openlayers/ol3/issues/3008

http://jsfiddle.net/mindsupport/gyumL4q6/3/

When i drag it or click, there is no property oData in object. Where is it ? But on  change listener work great.


texts: Clip raster with another raster (by extent) in Python

texts:
I'm trying to do raster algebra with two rasters, but when I transform the rasters to array I can't do the operation because the files do not have the same extent

example:

from osgeo import gdal
import numpy as np

IMG1 = gdal.Open('C:/IMG1.TIF')
IMG2 = gdal.Open('C:/IMG2.TIF')

img1 = IMG1.ReadAsArray()
img2 = IMG2.ReadAsArray()

img1+img2
>>>>> ValueError: operands could not be broadcast together with shapes (5353,3352) (5548,3232)

So.. if I see the gt of every img show that img1 and img2 have the same size of the pixel (100) but different extent

g1 = IMG1.GetGeoTransform()
gt1
>>>>> (482084.999999999, 100.0, 0.0, 5652715.00001374, 0.0, -100.0)
g2 = IMG2.GetGeoTransform()
gt2
>>>>> (475184.999999999, 100.0, 0.0, 5643715.00001535, 0.0, -100.0)

So.. I used this operation in order to get the intersect extent.. but I don't know how to apply

r1 = [gt1[0], gt1[3],gt1[0] + (gt1[1] * IMG1.RasterXSize), gt1[3] + (gt1[5] * IMG1.RasterYSize)]
r2 = [gt2[0], gt2[3],gt2[0] + (gt2[1] * IMG2.RasterXSize), gt2[3] + (gt2[5] * IMG2.RasterYSize)]
intersection = [max(r1[0], r2[0]), min(r1[1], r2[1]), min(r1[2], r2[2]), max(r1[3], r2[3])]


texts: Query to Dissolve Nearby Polygons on Hive

texts:
I am using the ESRI spatial framework for Hadoop, that extends Hive to use spatial types and operations.

My objective is to translate a set of simple queries on PostGIS into Hadoop, in order to reach horizontal scalability.

I have a grid with a count for each cell.

The objective of my query is to select all cells that have a count higher than a certain threshold, and group(merge) all cells that are together. For instance in this case, I would end up with something like this: 4 polygons.

To do this in PostGIS, I use a combination of ST_Dump and ST_SnapToGrid

CREATE TABLE exploded AS
SELECT
(ST_Dump(st_union)).geom
FROM  (SELECT ST_Union(ST_SnapToGrid(geom,0.0001))
 FROM grid where ptcnt > 'threshold) as q;

Unfortunately, none of these functions is available on ESRI's spatial framework.

I can perform the threshold filter, but I have no way of aggregating the nearby geometries based on the proximity (a trick perform by the grid):

create table exploded as select u as geom from (select geom as u from grid_cnt where ptcnt > 11467) as q;

Does anybody can think of a workaround (perhaps using Union)?


texts: Is there a way to map out PLSS (Public Land Survey System) TWP (township),RNG (range) and Section based off a field in google docs?  I was thinking if you had your first column with the PLSS area in it then it would be able to read the PLSS shapefile in ArcMap to adjust what section needs to be highlighted.  I want to be able to track my company's projects on map.  Right now they are all on a spreadsheet.  Any ideas to go about doing this?


texts: Can the ArcGIS coordinate change be wrong?

texts:
I am using my own Lambert Conformal Projection with following parameters.

False_Easting: 722.5,
False_Northing: 755.0,
Central_Meridian: -102.5,
Standard_Parallel_1: 20.0,
Standard_Parallel_2: 60.0,
Latitude_Of_Origin: 44.0,
User_Defined_Unit: 5000.0 m

I imported latitude-longitude data (WGS_1984) and converted its projection to my own projection above.

When I picked up one point having a coordinate of (lat=41.408,lon=-73.590)
, its coordinate is shown in ArcGIS as (xi=1167.6, yj=774.8).

However, when I calculate (xi,yj) of the point using WRF FORTRAN subroutine, it is (xi=1166.3, yj=774.6).

The error is about 1.3 grid (~6500 m). I think the difference between ArcGIS and WRF is too big.

Could you give me an advice for this problem?


texts: Raster terrain curvature analysis in QGIS

texts:
When QGIS performs a curvature analysis of the terrain raster, which curvature analysis does this represent i.e. the Gaussian curvature analysis, or Mean curvature analysis?

Understanding curvature filter of QGIS raster terrain analysis? is a similar question, where @whuber gave an extensive reply on how the curvature analysis is performed.

However what I still do not understand is whether or not the curvature analysis of the terrain raster is always a Gaussian curvature analysis, and never Mean curvature analysis.


texts: How to import OSM data into Postgres DB using osm2pgsql in Linux?

texts:
I'm using the following command to import an osm file into a Postgres DB using Linux:

osm2pgsql -v -a -U user -S default.style -d routing myfile.osm

I got the following error:

  Setting up table: planet_osm_point Problem reading geometry
  information for table planet_osm_point - does it exist? Error occured,
  cleaning up

I guess I have to manually create the schemas into db.

Can osm2pgsql be used to automatically create such schemas?
Otherwise, how can I create them?


texts: Difference between Trip and Route in Google Transit

texts: I am developing a transit webapp using Google Transit, while referring to documentation, i see two different datasets (trips.txt and routes.txt) but i can't clarify the difference between them.


texts: Clip raster with pyqgis and nodata value

texts:
I'm writing a script to clip many raster by a shape.
When a try the tool in pyqgis just for one raster, value around image are nadata, it's good !

Here is the code :

##test=group
##clip=name
##raster=raster
##shape=vector
##sortie=output raster

clip=processing.runalg('gdalogr:cliprasterbymasklayer', raster,shape,None,False,False,False,5,4,75.0,6.0,1.0,False,0,False,None,sortie)

After that I create a script with a loop on a folder. It's works but I have a little issue. Nodata value become 0 and I really don't know why.

Here is the new code :

##Fruition=group
##Clip_raster=name
##dossier_entree=folder
##dossier_sortie=folder
##parcellaire=vector

from qgis.utils import iface
import sys, processing, os, glob, numpy
from qgis.core import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *

fichiers = os.listdir(dossier_entree)
ndvibruts = glob.glob(dossier_entree+'/*.tif')
for ndvibrut in ndvibruts:

    Clip = processing.runalg('gdalogr:cliprasterbymasklayer', ndvibrut,parcellaire,None,False,False,False,5,4,75.0,6.0,1.0,False,0,False,None,None)

    sortie = QgsRasterLayer(Clip['OUTPUT'])

    extent = sortie.extent()

    width, height = sortie.width(), sortie.height()

    renderer = sortie.renderer()

    provider = sortie.dataProvider()

    crs = sortie.crs().toWkt()

    pipe = QgsRasterPipe()

    pipe.set(provider.clone())

    pipe.set(renderer.clone())

    nomndvibrut = ndvibrut[(ndvibrut.find('ndvi_')):(len(ndvibrut))]

    file_writer = QgsRasterFileWriter("{0}/{1}".format(dossier_sortie,nomndvibrut))

    file_writer.writeRaster(pipe,
                   width,
                   height,
                   extent,
                   sortie.crs())

Somebody can help me ?


texts: Creating polygon out of very large line feature class using ArcGIS for Desktop?

texts:
I am trying to create a polygon out of a line feature class in ArcMap 10.3.1.

I have tried using the advanced editing "construct polygons" tool and mesed with the clustering sizes but only get polygons in the very small tight corners of my would be polygon... I have also tried using the "feature class to polygon" tool but have gotten the same results.

I'm assuming at this point that my polygon in the main section is not completely closed and this is why it is happening the way it is. The problem is that I cannot check the topology due to the feature still being classified as a line (or I don't know how, if it is possible) other than visually, which I have tried to do for a few hours now to no avail.


texts: Problem with definition of categorical variables in dismo package

texts:
I'm working with a package dismo: Species Distribution Modeling and I have all my environmental variables as categorical. I used "factors" argument as follow:
mod <- dismo::maxent(x = pder,
                     p = pa,
                     path = maxent_out,
                     args = c(...),
                     factors = c("first var name", "second var name"...)
                    )

The problem is that from the html file that I get, all my variables are considered as continuous. At this point I don't know where is the problem.


texts: Programatically detect the format/version of an arc toolbox (e.g. 10.0,10.1 etc)

texts: Is there a way to tell programatically (from the command line) what format a given ArcMap toolbox was saved in - 10.0, 10.1 and so on?


texts: Why PyQt http request works in QGIS console, but not in a plugin?

texts:
I want to send http request in a QGIS plugin using PyQt.
I found in this post working syntaxes to do so, but something odd occurs.

The script given in the question works both in the console and in my plugin.
But the script given as an answer by German Carrillo, technically better, works when executed in QGIS python console, but doesn't as soon as it is integrated in a plugin.

So this works in the console, but doesn't when integrated in a plugin script :

from PyQt4.QtCore import QUrl
from PyQt4.QtNetwork import QNetworkAccessManager, QNetworkRequest

def request_download(manager):
    url = QUrl("https://www.facebook.com")
    request = QNetworkRequest(url)
    manager.get( request )

def handle_download(reply):
    print reply.readAll()

manager = QNetworkAccessManager()
manager.finished.connect( handle_download )
request_download(manager)

Considering the logs I have put in my script to locate the issue, it appears the handle_dowload function is never called. So it means the QNetworkAccessManager doesn't send his finished signal.

Is there known reasons why something would work in the console but not in a plugin ?


texts: How to filter in the attribute table within an attribute name?

texts: I have to create a map with the corridor information in it (in German: Flurstücke). For that, I have one shapefile with the lines and another shapefile (points) with the information about the corridor numbers (and a lot more information). Now I would like to filter the the point-shapefile and delete all the points that don't have the corridor number or land use. Is it possible to filter the attribute table within an attribute name (for example within the attribute name Layer)? To be clear, I want to delete everything that doesn’t have information about the corridor number or the land use. How would be the formula for that?


texts: Drawing a straight line and finding neighbouring area

texts:
I have a shapefile showing district boundaries, and the coordinates of cities. I need to:

Draw a line between two coordinates e.g. draw a line between the centre of London city and the centre of Manchester city.
Select the districts where the line intersects
Select the neighbouring districts (districts directly beside the districts in 2.)

How can I do this?

I am sorry I do not have any figures to show.


texts: GeoJSON won't show in IE, but works fine in FF

texts:
I’m not sure if this question isn’t too specific for Stack Exchange, but I really can’t find a solution, and even extensive web search hasn’t brought up any help. There are similar, but different questions here on Stack Exchange, but they also didn't help me.

Basically, I am trying to display a simple, rather small geoJSON file in an OpenLayers project, just about 20 polygons (a small showcase, the focus is more on tiled imagery than on the geoJSON). I created the geoJSON with QGIS from an existing shapefile and cut off all attribute data (I also validated it online to be sure)

This is the actual website (containing the geoJSON showing city districts, and two tiled planning documents that are disabled by default):

http://www.anderka-stadtentwicklung.de/maps/penzberg/webmap.html

(sorry for the abysmal performance, my webhost is moving his server farm this week)

A quick link to the actual JS map code:

http://www.anderka-stadtentwicklung.de/maps/penzberg/webmap_code.js

The polygons in questions are shown in orange, at least in Firefox. But soon I noticed that the geoJSON wasn’t visible when I viewed the map in IE9 (neither in normal mode, nor in any of the various legacy and development modes). There is no error message, and everything else works fine. I already checked all the usual suspects (trailing commas, projection, proxy problems etc.), none of it helped. JS isn't really my strongest skill, but I think my code is more or less valid, because Firebug won't point out any glaring errors, and the map is showing up correctly in FF...

The JS code itself mostly stems from a former project of mine (loosely based on an example from Erik Hazzards great book), plus various snippets I gathered up here on this page, some of the Openlayers tutorials and different Q&A sites on the web.

Right now I am also using the OpenLayers 2.12 full feature .js file to avoid any unintended side effects due to missing features (the long-term goal is to create a slim custom openlayers.js, of course).

So my question is, in short, why won’t the geoJSON show in IE? What am I doing wrong here?

I would be deeply grateful if someone could at least provide me with a hint where to look for a solution… at the moment I’m completely lost.


texts:
I have daily SO2 data (the worldwide data) and after some processing I have to sum 30 rasters.

What the tool can you recommend to use for this goal within ArcMap tools?

I tried Raster calculator ("rast 1"+"rast 2"+...), but it doesn't give me the raster with the sum of pixel values.

Also I tried Image analysis (Mosaic tool), but still no the result I need to obtain.

Furthermore the summarizing by Raster calculator causes the missing of some values in different parts of the map.


texts: Converting GRIB to Geotiff with gdal_translate in Python

texts:
I want to extract each band from a GRIB (gridded binary from http://data-portal.ecmwf.int/data/d/interim_daily/) in a single Geotiff image. I´m using Python with gdal_translate in a QGIS backround, but the dgal_translate command doesn´t work.

The import of the GRIB-data works well, but in my loop about all bands, where I want to execute the gdal_translate command I get an error.
The loop also works well, but it tells me, that the algorithm for the GDAL command is unknown (I imported the tool in my script).

I get the information (in the QGIS console) while running: "Algorithm not found" (I get this information for each band, but the program continues).
There is no other explicit error warning.
My code:

import os
from osgeo import gdal, ogr, osr from

filename='C:\\Users\\mypath\\test1.grib' rasterInput =
gdal.Open(filename, GA_ReadOnly)

if rasterInput is None:
    print 'Grib import failed'

print 'Grib imported'

bandCount = rasterInput.RasterCount print bandCount

format = "GTiff" driver = gdal.GetDriverByName(format)

for a in range(bandCount):

    outputFile = 'C:/Users/mypath/Geotiffs/%s.tif' %(a+1)

    translate = processing.runalg('translate, -b %s, -co COMPRESS=NONE, -co BIGTIFF=IF_NEEDED, -of GTiff, %s, %s' %(a+1, rasterInput, outputFile))

del bandCount

I copied and pasted this command from the QGIS tool raster --> conversion --> translate, but it just works for one band not for all in range automatically.

How can I do this per command-line?

I´m a beginner at Python and GDAL and I can´t find any satisfying result.

I hope that is explicit enough for some answers.


texts:
I have been assigned a task to create a THREDDS server so our lab can start uploading data online. I managed to learn how to download and work with Tomcat/THREDDS and established our own local server. I have also been feeding data into the server without any problems so far, but I am confused on how can I make my catalog available for the public.
The local server can only be accessed from our computer or through remote access as far as I understand; but how does one make it available for open access so users can access and download data? I know databases such as Hycom have a THREDDS catalog page in their domain.
How can I make this possible?


texts: Changing data from ETRS89 30N (Spanish local CRS) to WGS84 to use GPS tracker in Spain

texts:
I have made up a number of maps in QGIS in order to complete a small mapping project in North West Spain. The data was downloaded from IGN:

The vector files downloaded were from: centrodedescargas.cnig.es/CentroDescargas/index.jsp , from the 1:25,000 vector maps in .dgn format.

These are referenced to CRS system ETRS89 30N. I've attached an example of one of the maps.

I have a few questions (I'm not very experienced in GIS so please bear with me...)

How can I change the grid system on my map to match with WGS84 when my layers are in ED50 30N and actual grid is in ETRS89 30N, ensuring everything is properly calibrated and matched up?
Can I use British grid on my Garmin Etrex 10 (handheld GPS) in Spain, considering they're both in ETRS89?


texts: How to create an ArcGIS toolbox programmatically?

texts:
Is it possible to programmatically create an ArcGIS toolbox?

The aim is to write a plugin that works in ArcGIS and have a single source for the tool parameters which is then compiled to an Arc toolbox app.

Alternatively is there a good workaround for this?  I note data interoperability in Arc fires up a second dialog from the toolbox dialog: can this be done from the tool validator?  As such a dialog could obviously be compiled by an external script, this (almost) fulfils the same purpose.

CLARIFICATION this question is about Arc toolboxes only, as I regard programmatic creation of QGIS/command line interfaces to be a solved problem already.


texts:
I'm currently trying to use QField on Trimble TDC600s. When I try and enable GPS position display on the map, the display clears off and leaves a single dot on it.  Every other feature is working as intended including the display of position information at the bottom of the screen.

Please note: I have tried reinstalling as that seemed to fix this for at least some others in the past.  This also occurs on multiple units.  Additionally when I long press to turn on and off positioning, I get an error message that QField does not have permissions to positioning when turning it off, but not on.

QField version 1.3.1
QGIS version 3.10. Data layers are in a geopackage on the Trimble.


texts: Making line from set of GPS points using ArcGIS for Desktop?

texts:
I need to create the trajectory of respondents according to the their collected GPS data (including X and Y).

I am familiar with python programming, but I am a beginner in ArcGIS.

I want to write a script to automate this process.


texts:
I'm asking this question for a student I have that can't open shapefiles in QGIS.

Whenever she tries to open shapefile in QGIS 3 she get this error message:

When I checked her library where the shapefile stored, I saw this:

these files, that on my computer are shapefile, on her computer are AutoCad.

This is how it looks in my computer:

My question is, how can I change it so she can load shapefiles into QGIS?


texts:
I am trying to create a vector grid (polygons) but the tool Vector>Research Tools>Vector grid lets QGIS crash all the time. I tried different grid extents (same as canvas, by another layer), resolutions, opened a blank project ... but it keeps on crashing.

I am using the newest QGIS 2.18.3


texts: Using Define Projection and Projection tools in ArcMap

texts:
When I uploaded my maxent output file in ArcMap (it's in ASCII format, a distribution map of a country) it showed me the the spatial reference is unknown. I will create binary maps later (suitable, unsuitable habitat) and calculate the area. So I have to make the file ultimately in a UTM coordinate since its unit in meters. I'm confused which approach of the following is accurate, shall I use Define Projection tool to define it to WGS84 (DD unit) first -which is the original coordinate system of the data I entered into maxent- then use Project tool to change to a projected coordinate system (UTM, meters unit)? Or I can immediately from the define tool choose a UTM coordinate?
If the first approach is the accurate one. I did tried it but when I reach the project tool step it doesn't recognize the file in the input. I tried exporting the file to a tiff the same issue. When I saved the file as a layer (*lyr) it did recognize it but showed me this: ERROR 000840
The value is not a Feature Layer. ERROR 000840: The value is not a Feature Dataset.


texts: Process PyQGIS 3.8 qgis:zonalstatistics and save output to CSV without geometry

texts:
processing.algorithmHelp('qgis:zonalstatistics') prints out:

----------------
Outputs
----------------

INPUT_VECTOR:  <QgsProcessingOutputVectorLayer>
        Zonal statistics

but I'd like to redirect the output to a CSV file with only the id column of the INPUT_VECTOR and the output of the Zonal statistics algorithm. Is this possible??

Is there some way to redirect the output?


texts: Why is the display of QgsMessageBar delayed?

texts:
When a function is called - a file is tested for existence.  This checksum pushes a messageBar() object for both circumstances (True/False), yet the actual bar does not
display in MapCanvas until after the subsequent exec_() is called and finished running.

def callUtherPendragonsSon(self):
    holyGrailFile = 'C:\ALongWayFrom\camelot.csv'
    if os.path.exists(holyGrailFile):
        self.iface.messageBar().pushMessage("Notice", "You will need many coconuts...",\
                                            level=QgsMessageBar.INFO, duration=3)
        beginTheJourney = gatherProvisionsAndMen()
        beginTheJourney.exec_()
    else:
        self.iface.messageBar().pushMessage("Warning", "There is no grail :( ",\
                                            level=QgsMessageBar.CRITICAL, duration=3)
        return

Is there a means to assure the message bar is pushed prior to calling the exec_() ?

UPDATE (1):

As it turns out, a command in the external Ui gatherProvisionsAndMen()'s init  is causing the trouble.  It is a data-thirsty function call, and seems to supersede the messagebar (even though it is executed thereafter).   The messageBar itself is intended to preliminarily notify the user of this nested 'lengthy process' - but it seems such a notice needs to take place in another manner...

Perhaps I hold a misunderstanding in the logistics of programming but,

wouldn't the call to push a QgsMessageBar() need to complete before moving forward and executing the external Ui?


texts: pyproj EPSG:4326 converts degrees to radians

texts:
Case

My input is a point given as x, y and an SRS name which specifies how to interpret the numbers.  My goal is to convert these into lon, lat in degrees.

I convert this by using this piece of code:

projection = pyproj.Proj(srs_name)
lon, lat = projection(x, y, inverse=True)

A typical input is coordinates in meters and the name of an SRS which has meters as unit.  The code above works fine in this case.

Problem

Sometimes the input is the SRS name "EPSG:4326", and the input coordinates then are given in degrees.  From what I understand about EPSG:4326 this is correct, too.  Since the input already is in degrees, the conversion done by pyproj should not change anything, i. e. projection(x, y, inverse=True) should return a lon, lat with lon == x and lat == y.

Unfortunately this is not the case.

I figured by now that the projection = pyproj.Proj('EPSG:4326') converts the input from degrees to radians.  So if I use the inverse=True variant, it converts radians to degrees.  Since my input is not given in radians, the result will be nonsense (in degrees).

What I Tried

I played around with the kwargs radians=False and preserve_units=True.  I tried giving it to the pyproj.Proj() and to the call of the conversion (projection(x, y, invers=True, radians=False)).  None of the combinations I tried had any influence.  The inverse=True conversion always assumes radians as input.

I also tried figuring out if there is some kind of global switch in pyproj to make it use degrees instead of radians.  I didn't find any.

I'm using pyproj 2.2.0.  From the git log of the project (as in Github) I understand that there have been some changes made in the recent past concerning the kwargs radians.  Version 2.2.0 could be older than these, and the changes could be a fix for my case.

I also found on https://epsg.io/4326 that the unit of EPSG:4326 is specified as UNIT["degree",0.0174532925199433, ...] and this number is the factor between radians and degrees.  This now looks like the unit of this EPSG in fact is radians.  Should this be the case, then I wonder why my input is in degrees instead of radians if the SRS is EPSG:4326.  Maybe the client (providing the input) does it wrong?

Question

How do I have to change my use of pyproj in order to not run into the problem anymore?