texts: Determining longest line segment within polygon that passes through its centroid?

texts:
A PhD student approached me recently asking how one would determine the longest straight line that passes through the centroid of a polygon, the output being a polyline rather than just a table of numbers. The polyline would be within the polygon starting at the edge at one end, passing through the centroid, and then end at the opposite edge.

Does anyone know how to calculate this?

I'm surprised (but I guess I'm showing my ignorance) that this is not some interface in ArcObjects.

The polygon represents the crater edge of a volcano with no islands so the polygon can have an irregular shape.


texts: Is there any way to view maps in "perspective" with ArcGIS (ArcMap do not)?

texts: I need to work with ArcGis by the amount of data that I use and the answers I need, but would like to see my maps in perspective and I can't find a way. Does anyone knows how to do it?


texts: Errors when running python gdal.grid with basic parameters

texts:
There are a few problems when I run the following code.

Image comes out flipped vertically, even though lat lon is input correctly
I have to run the code twice for it to run properly (so wierrd)
I also get the following errors:

    ERROR 4: `result.tif' not recognized as a supported file format.
    ERROR 4: Unable to open result.tif to obtain file list.
    ERROR 1: Deleting result.tif failed:
    Permission denied

I use the following data in the points.csv:
    lat lon value
    39.545  -75.78  0.001034494
    39.555  -75.78  0.001033509
    39.565  -75.78  0.001032392
    39.575  -75.78  0.001001829
    39.585  -75.78  0.00094428
    39.595  -75.78  0.000940016
    39.605  -75.78  0.000897172
    39.615  -75.78  0.000896742
    39.625  -75.78  0.000896021

And here is the code...
# create the csv and vrt file
vrt_file = """<OGRVRTDataSource>
    <OGRVRTLayer name="points">
        <SrcDataSource>points.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <GeometryField encoding="PointFromColumns" x="lon" y="lat" z="value"/>
    </OGRVRTLayer>
</OGRVRTDataSource>"""
with open("points.vrt", "w") as f:
    f.write(vrt_file.strip())

# perform gridding and RGB rastering
# for some reason this has to run twice????
gdal.AllRegister()
gdal.Grid('result.tif',
               'points.vrt',
               format='GTiff',
               outputType=gdal.GDT_Float32,
               algorithm='linear:radius=0.5',
               zfield='value')

# plot, comes out flipped on accident
im = gdal.Open('result.tif').ReadAsArray()
im = plt.imshow(im)


texts: Including NoData zones in zonal statistics table

texts:
I am running the zonal statistic tool on a raster with a shapefile that consists of several regional polygons. I would like the zonal statistic to include the regions that are completely NoData in the final table instead of just completely removing it. Is there a way to do this or am I limited to the Data or NoData options that are available?

My current code is simply:

arcsa.ZonalStatisticsAsTable(carPath, 'CARUID', os.path.join(tmpDir, "masked.tif"), os.path.join(tmpDir, "table"), 'DATA', 'ALL')


texts: Is it possible to find if a value is present in Array without looping in mapbasic?

texts:
How to find a  value is present in array without looping ?

e.g.
For i=1 TO ubount(Array)
  Find if value is present
Next

or some more direct statement like

Array.contains(somevalue) like in java


texts:
I am trying to georeference a raster map from Nepal.
What CRS should I use? PROJ.4 or WKT form.
I will use the UTM grid to base my reference points on.

On the map there is reference to:
Scale 1:50000
Everest 1830 Spheroid
Modified UTM Projection
Origin Long 87E Lat 0N
False Coordinates at Origin Easting 500000 Northing 0
Scale Factor at Central Meridian 0.9999

Centre of map is about 460000E 3080000N


texts: How to process multiple images in Grass GIS

texts:
first of all, I have to mention that I'm almost new working with Grass and I'm using Grass 6.4.3 on Linuxmint 16. I need to process a large number (~300) of Landsat images. I have downloaded the images from the USGS, I have all of them unzipped, and one folder for each scene containing one file for each band and the metadata file. I need to read/import all the bands of each scene, apply them the radiometric conversion and the atmospheric correction (i.landsat.toar), then do some calculations (r.mapcalc) and finally write a new raster file with the results of these calculations (r.out.gdal?). I can do it manually, importing all bands for each image and doing all the process one by one. That could take me a long time. Much better would be to run a batch job. Could you help me to do it?
I've been searching and reading a lot on the web, but I didn't find a solution.
So, my problem is: how can I do to read one file, process it, save it, and then read the following file, etc. I have a background programming (R and other languages, not in Python) but I couldn't understand how to perform a batch job inside Grass. Is there a step-by-step tutorial explaining how to do it? Is it possible to read external files and folders in Grass? or I have to import all the raster files to the mapset first and then run the batch job?

Thanks a lot!

EDIT: I think I have reached to something useful, now I have this code:

#!/bin/bash
# Directory where the images are stored
ImgDir="/home/lbianchi/Images/"

for Dir in $(find $ImgDir* -maxdepth 0 -type d ); do
cd $Dir

# Import all files(bands) in de directory
for files in $(find . -iname '*.tif'); do
  thefile=$(basename $files)
  importedname=${thefile:0:24}
  fullpath=$PWD"/"$(basename $files)
  r.in.gdal input=$fullpath output=$importedname
done

# once the bands are imported, apply i.landsat.toar
inputname=${thefile:0:23}
correctedname=${thefile:0:22}"toar_DOS4_"
metadata=${fullpath:0:97}"MTL.txt"
i.landsat.toar input_prefix=$inputname output_prefix=$correctedname metfile=$metadata sensor=tm5 method=dos4
echo $inputname"--"$correctedname
cd ..
done

echo "the end"
exit 0

It does almost what I need, can you help me to improve it?
And more important, is there a way to do the same process but without importing ALL files to the Mapset. Otherwise I'm thinking in deleting the imported files once the process is finished.
Thanks again.


texts:
How can I edit selected Attribute only with python 2.7

I want to select points in arcgis and then when I click on the button created by my python script I want the value of the selected features change automatically.

when i run the code, all the features with the value 99 changes to 4 and I want to only change the value of the features that are selected in my current mxd in ArcGIS Desktop

this is my script

 import arcpy import pythonaddins

 class EditFild(object):
     """Implementation for Tools_addin.EditFild_1 (Button)"""
     def __init__(self):
         self.enabled = True
         self.checked = False
         print('int')
     def onClick(self):
         features = arcpy.UpdateCursor(r"D:\\SHP_test\\Point.shp")
         for feature in features:
             if feature.x == 99:
                feature.x = 4
                features.updateRow(feature)
         # Works fine
         print("done")
         del feature,features


texts: Calculating flowpaths downstream of a point

texts:
I am assessing the flood risk posed by pipeline failure.  I need to calculate the likely flow-paths for flood water downstream of critical locations along the pipeline using GIS.

I'm familiar with ArcGIS and QGIS/GRASS/SAGA methods for calculating channel networks and catchments upstream from a point but what would be the best method for assessing downstream flow-paths originating at a point location?

Ideally any method would be fairly automated as a number of locations are being assessed.

I have ArcGIS 9.3 and QGIS (latest version) available to me.


texts: How to join points into overlapping polygons layer?

texts:
I'm using QGIS and I am relatively new to it. I have a points layer. Each point is the location of a conflict event and the attributes are the type of conflict, year of the event, number of deaths etc. I also have a polygon layer. The polygons are the areas inhabited by different groups which can overlay. Two final results work for me:

I want to have the number of events BY TYPE, YEAR, etc for each polygon.
It is also useful if I could have a table with the points and a column with the polygon attributes. In this case, if a point falls into two different polygons, I will need that point to appear twice once for each each polygon.

Please let me know whether there is a way to do at least on them in QGIS.


texts:
I am trying to import an Excel table with around 10,000 XY points into ArcMap 10.2. I am using State Plane NAD 1983 Alabama West (US Feet) as the Coordinate System in my map data frame. I then add the excel table to my Table of Contents in the map then select Display XY Data. Next, I choose the same coordinate system as my data frame then select OK. Once it creates the Events layer for my points they are not in the correct location.

Should I switch everything in my map to WGS84 or is there a solution for my points to be added onto my current map with the state plane coordinate system?


texts: Configuration problems with GeoServer 2.6 under Ubuntu

texts:
I installed GeoServer on Ubuntu Linux 14.04 using the OS-independent binary that comes with Jetty. I am running it using the bin/startup.sh file. I am trying to create an ImageMosaic. When I try, I get this error message in the browser console:

ReferenceError: Wicket is not defined

I also get odd error messages in the console here and there, complaining about not finding JQuery - and yet, according to Firebug, the JQuery JS and CSS files get found and download OK.


texts:
Can u help me in creating a trigger and function so that every time i change the flag in table A the buffer distance should get changed and also its geometry field also , ie the polygon area?

I have two tables named A and B.

Table A has the fields GID (primary key), FLAG, Buffer_distance, and the geometry field.
Table B has the fields GID and buffer_distance, and has only three records as shown below:

GID(PK)   Buffer_distance(Double Precision)
1     1000
2     2000
3     3000

Table A's structure:

gid(PK)  buffer_distance(Double Precision)  flag(Text)  the_geom( Geometry)
1           200              1          the_geom
2            100             2          the_geom
3            100             3          the_geom
4           500              3          the_geom
5           300              2          the_geom
6           899              1          the_geom

What I want is the value of flag in table A to be updated from table B, and as soon as I updated the value of flag in table A the buffer_distance in A should also get replaced along with the geometry.

If what you want is to update buffer_distance in Table A based on the flag value (ie, buffer_distance=2000 when flag=2),  can anyone  write me an SQL case statement that will accomplish this.
along with this the_geom Filed the geometry should also change with the new buffer_distance.
i used the following code.
help me to correct the same

create or replace function update_point_buffer() returns trigger as
$$
 begin
 -- delete
IF (TG_OP = 'UPDATE') THEN
UPDATE tablea
    SET buffer_distance = (SELECT tableb.buffer_distance
                        FROM tableb
                       WHERE
                       "tableb.gid" = "tablea.flag");

   END IF;
   return NEW;
END;
$$ LANGUAGE plpgsql;

CREATE TRIGGER trg_upgradeBufferTrigger
        AFTER UPDATE ON tablea
        FOR EACH ROW
        EXECUTE PROCEDURE update_point_buffer();

i have used this query bit it has given this error ERROR: syntax error at or near "v_output" LINE 2: v_output integer;

create or replace function update_point_buffer(v_pointflag) returns number as---------------------Return buffer distance v_output number ; begin

SELECT buffer_distance into v_output FROM point_buffer WHERE gid = v_pointflag;

return v_output; END;

CREATE TRIGGER trg_upgradeBufferTrigger-----------------------------------Update table from B to A AFTER UPDATE ON point FOR EACH ROW

begin

update point_buffer set buffer_distance= (select buffer_distance from table B where gis= :new.gis)

where gis= :new.gis

end;

can any body correct it for use.


texts: Placement of a fixed-radius circle on a sphere to maximize point containment

texts:
There are many papers on the "Circle Placement Problem" in 2D Euclidean space, but I haven't been able to find any papers discussing this problem on the sphere. The problem is: given a set of points on a sphere and a fixed radius for a circle, find the location to place the circle that will cover the maximum number of points.

An example of a paper that discusses this problem in 2D is "On a Circle Placement Problem" by Chazelle and Lee. As discussed in this paper, this problem could be used to find the best location for a radio transmission tower with a fixed range given the locations of nearby cities.

Any thoughts you may have would be welcome. Thanks!


texts:
I have a raster (.sid extension) file. I want to digitize roads automatically in ArcMap, QGIS, etc. Digitizing by hand would be too time consuming.

How can I do this ?


texts: How to get features of another point layer when geomequals ? (Using Default Values in Attributes Form Qgis 3.2.1)

texts:
I have 2 different point layers with attributes fields as shown in the picture below. They both have a different amount of features. Now i want that all features of layer B gets filled with some feature attributes of layer A as long they are on the map at the same position. (I used the snapping function to be sure that some positons between the 2 Layers should be equal)

So was testing to use default values for each field in layer B and failed. For example:

Layer B Field B1 - Default value: Case when geomequals(A,A1) then 1 end

The output I got here on all features of Layer B for fields B1 - B3 using that expression was always just NULL instead of 1. (The 1 was just for testing)

Anyone here got an solution for this issue and an additional expression for 1 which refers to the attribues in layer A for the fields A1, A2, A3 ?


texts: How to handle Molusce ValueError: shape mismatch: objects cannot be broadcast to a single shape?

texts: I have a problem in QGIS Molusce plugin. I opened two raster files represeting LULC of 2 years, and roads as a spatial variable. The result of geometry check was positive. On Area changes tab clicking Update tables results in a Python failure message: "ValueError: shape mismatch: objects cannot be broadcast to a single shape". Does anyone have an idea what causes this failure message? (I tried this in different versions of QGIS - 1.8, 2.2, 2.6 resulting in the same message.)


texts: Position of a point changes by zoom in QGIS ? (OpenLayers plug-in doesn't show same position on different zoom levels ?)

texts:
I have a point layer in QGIS. I was just inspecting the points and I've realized points position changes too much as I zoom in or zoom out. Does anyone know why this is happening?

I'm sure Google layers plug-in is working well because many people using it everyday. Also I've checked my point layer is referenced to true coordinate system.

My question is quite straight forward and maybe stupid I don't know but I couldn't figure out why so I decided to ask it here.


texts: Why does updating `rgdal` package in R change output of `sp::CRS` (+towgs84=0,0,0)?

texts:
I am not a GIS expert, and this change caused quite some errors in our analysis. So my question is, why this was changed: From rgdal 1.4.8 to 1.5.16, the +towgs84=0,0,0 disappeard in the output of sp::CRS.
I hope for an answer, which is comprehensive to a non-GIS-expert.
Reproducible example if you have NO sp and NO rgdal package installed:
packageurl <- "http://cran.r-project.org/src/contrib/Archive/sp/sp_1.4-1.tar.gz"
install.packages(packageurl, repos=NULL, type="source")

packageurl <- "http://cran.r-project.org/src/contrib/Archive/rgdal/rgdal_1.4-8.tar.gz"
install.packages(packageurl, repos=NULL, type="source")

packageVersion("sp")
# [1] ‘1.4.1’
packageVersion("rgdal")
# [1] ‘1.4.8’

sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
# CRS arguments:
#   +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0

So everything works, that is in our case
spdf <- sp::SpatialPointsDataFrame(
  coords = mission_coords[, c(2,1)],
  data = mission_coords,
  proj4string = sp::CRS(
    "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
df_in <- sp::over(spdf , departement_boundaries , fn = NULL)

are ok. (I have problem to make a small reproducible example - the departement_boundaries from Zurich are open source, but its a huge string from dput.)
Now, update both packages
install.packages("rgdal")
# also installing the dependency ‘sp’
packageVersion("sp")
# [1] ‘1.4.2’
packageVersion("rgdal")
# [1] ‘1.5.16’

and in any case +towgs84=0,0,0  is missing now:
sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs",
        doCheckCRSArgs = TRUE)
# CRS arguments: +proj=longlat +datum=WGS84 +no_defs
sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs",
        doCheckCRSArgs = FALSE) # ... +ellps=WGS84 is added...(?)
# CRS arguments: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

Any ideas for an explanation to a GIS-newbie?

@EdzerPebesma’s answer helped me find R spatial follows GDAL and PROJ development to get information:
x <- sp::CRS("+init=epsg:4326")
cat(comment(x), "\n")

One comment from that site is really unclear to me (perhaps how this relates to +towgs84=?):

rgdal generates many warnings about discarded PROJ string keys,
-> See e.g. Warnung in showSRID(uprojargs, format = "PROJ", multiline = "NO")
Discarded datum WGS_1984 in CRS definition, but +towgs84= values preserved
from here?
intended to alert package maintainers and script authors to the need to review code.
-> To me it is not clear, what I have to do, see the link above.
It is particularly egregious to assign to the CRS object projargs slot directly,
and this is unfortunately seem in much code in packages.
-> Again, what should I do instead?


texts:
A viewer is being developed with the Api Javascript of Esri 3.21 and several queries are made through QueryTask to some services that we have published in ArcGIS Server. It has been observed that the service response is somewhat slow.

For example, to perform this query on my machine (Window 8.1, Intel(R) Core(TM) i5-5200U CPU @2.20GHz) Google Chrome(Versión 61.0.3163.100):

https://186.3.64.173:6443/arcgis/rest/services/ServiciosDePrueba/RecaudacionDeuda/MapServer/2/query?where=PersonasEm%3E0&text=&objectIds=&time=&geometry=&geometryType=esriGeometryEnvelope&inSR=&spatialRel=esriSpatialRelIntersects&relationParam=&outFields=clavepre%2C+emitido%2C+pagado&returnGeometry=true&returnTrueCurves=false&maxAllowableOffset=&geometryPrecision=&outSR=&returnIdsOnly=false&returnCountOnly=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&returnDistinctValues=false&resultOffset=&resultRecordCount=&f=html

The query takes approximately 30 seconds and the records returned are 85,100 and the total service entities are less than 101,000.

We do not know if you can solve the problem in ArcGIS Server configuration to optimize this type of queries.


texts:
I'm a basic/novice user with GeoServer. I've installed it and have loaded a sample shapefile, which I can view using the Layer Preview option.

Now I want to add this shapefile (as a WMS layer) to a simple OpenLayers map using the methodology listed at http://dev.openlayers.org/docs/files/OpenLayers/Layer/WMS-js.html:

var wms = new OpenLayers.Layer.WMS("NASA Global Mosaic",
    "http://wms.jpl.nasa.gov/wms.cgi",
    {
    layers: "modis,global_mosaic",
    transparent: true
    }, {
    opacity: 0.5,
    singleTile: true
});

How do I find the URL which corresponds to one of the layers in the screenshot above?

If I use the Preview Layer option and choose JPG as the format, the URL is in the form:

http://localhost:8080/geoserver/cite/wms?service=WMS&version=1.1.0&request=GetMap&layers=cite:States2006 ... rather than the /wms.cgi form listed above.

Thanks for any advice on getting started with GeoServer and OpenLayers.


texts: Identifying row and column location of a raster image using ArcGIS 10.1

texts:
The solution already available online gives a python syntax error

  ERROR 000989

  Python syntax error: Parsing error SyntaxError: invalid syntax (line
  1)

when I write

$$rowmap + 1

in map algebra expression.


texts: Is there a way to give the drawing red polygon a higher alpha channel

texts:
Is it possible to set the alpha channel a bit higher when drawing? I don't mean the alphachannel of a already finished polygon, that's possible for sure!

I mean the red area. It's really stressful because I need to see hwat is beyond that polygon and in many cases my painting is harmed.

So is there any way to do so?


texts: How to switch to edit mode vector layers, in pyqt application based on qgis?

texts: I do QGIS based application, I need to implement the ability to modify the geometry of vector layers in manual mode (editing nodes), how to do it?


texts: Passing additional arguments to ST_MapAlgebra callback?

texts:
I've written a custom callback function for ST_MapAlgebra and I want to pass additional arguments to it. The function is pretty useless - it just generates a raster filled with normal random values according to a user-specified mean and standard deviation.

I tried including the extra arguments in the string, like so:

...
ST_MapAlgebra(
    b1.rast, -- raster to operate on
    1,
    'generate_random_raster(double precision[][][], integer[][],  ''3'', ''1'')'::regprocedure
  ),
...

And like this:

...
ST_MapAlgebra(
    b1.rast, -- raster to operate on
    1,
    'generate_random_raster(double precision[][][], integer[][],  VARIADIC ARRAY[3,5]::text')'::regprocedure
  ),
...

It works fine if I hard code the values I want to use, so the callback itself isn't a problem.

How do I pass additional arguments to a custom ST_MapAlgebra function?


texts: Migrating from PostGIS to Oracle Spatial

texts: We are now using a solution based on PostGIS extension of Postgresql database for various spatial queries. But we are now thinking about data security. This is why we have a plan to migrate existing PostGIS database to Oracle spatial. I do not know whether it is good idea or not. Can anyone suggest me what will be the best for us and what kind of challenges we have to face for the migration and queries?


texts: Connect points in two layers using lines in QGIS 3.x

texts: I'm using QGIS 3.2.1 and I have two Layers containing point objects, Is there any plugin for this version of qgis that I can use to select two points and then connect them using a line and save the line in a separate layer? I think I'm searching for something like ReconstructLine, But apparently it's not compatible with this version of qgis.


texts: Calculating elevation profile of a trail?

texts: What would be the best way to calculate an elevation profile of a trail remotely (without using altimeter/gps tracking your steps)? For example, to do it in the most rudimentary way you might take a topo map and measure distances and elevation changes along the course of a trail. Can you do this with USGS topo maps or do you need a DEM? What would be the best data source? I don't have access to Arc, just QGIS.


texts: GRASS GIS: r.out.gdal raises an exception and hangs if it's called multiple consecutive times

texts:
I'm using a small python routine (grass.script, pygrass) to handle several operations on sentinel2 imagery. I'm calling r.out.gdal inside the script and a very annoying problem occurs: when it is called multiple times (performing the exportation of several rasters) it raises an exception and hangs.
I've done several tests to check if the problem lies in my implementation, but the results tell the contrary: the hanging happens randomly at different points of the execution for the different tests (for some tests it manages to write all the maps but hangs instead of proceeding to the following operations, for some other runs of the same test script it only managed to write 2 or 3 maps and hanged). I then tried to use python time.sleep() and the results got better: it managed to conclude the operations for several runs, but sill hanged for 2/3 runs (using the same small testing script).
I'm pasting the text of the exception below.
Does anybody know why it happens and how to deal with the issue?

r.out.gdal complete. File
</Volumes/Maxtor/GISDBASE_histogram_analysis/clipped_maps_sectors/sentinel/20190223/rhow_20190223_B01_southern_wetlands_sector.tif>
created. Checking GDAL data type and nodata value... Exception in
thread Thread-1 : Traceback (most recent call last):   File
"/Applications/GRASS-7.8.app/Contents/Resources/lib/p
ython3.8/threading.py", line 932, in _bootstrap_inner
self.run()   File
"/Applications/GRASS-7.8.app/Contents/Resources/gui/w
xpython/core/gconsole.py", line 162, in run
self.resultQ.put((requestId, self.requestCmd.run()))   File
"/Applications/GRASS-7.8.app/Contents/Resources/gui/w
xpython/core/gcmd.py", line 603, in run
self._redirect_stream()   File
"/Applications/GRASS-7.8.app/Contents/Resources/gui/w
xpython/core/gcmd.py", line 632, in _redirect_stream
self.stderr.write(line)   File
"/Applications/GRASS-7.8.app/Contents/Resources/gui/w
xpython/core/gconsole.py", line 299, in write
value = int(line.rsplit(':', 1)[1].strip()) IndexError : list index
out of range


texts: (Auto)Convertion from point-coordinates (postgres)

texts:
When loading data from my postgres-db (where I can only read from, so not my own db) the field with the point-coordinates like ..
'POINT(390002.1610814564 6001433.406049595)'

.. are getting (auto)converted into ..
0101000020E964000022936245FEDA1741A4A29D5E0DE85641

Working with coordinate-systems is rather new to me, so can anyone tell me how this converted coordinate is called or how I can convert it back into the previous format?
Don't know if this information is helpful, but I'm using django (python) with the following code:
try:
    with connections['db_second'].cursor() as cursor:
        cursor.execute("SELECT * from ...")
        row = cursor.fetchall()
except Exception as e:
    cursor.close()


texts: Is ArcSDE no longer needed at ArcGIS 10.1?

texts:
Through my frustration of not finding the post-installation wizard with ArcSDE 10.1 and being unsuccessful with the command line tool to start an SDE service on our SDE databases, I found I was able to access the databases and copy/past/add features/delete features/etc. all through direct connections to our SQL Server 2012 machine. This was WITHOUT ArcSDE running.

So... The question is, why do I need to go through the trouble of trying to get the SDE services to run if I can manage our databases through the Desktop tools? I was even able to create an enterprise geodatabase and verified that all of the necessary business tables were created.
I have even recreated all our map services in ArcGIS Server 10.1 with direct connect and the web editing is perfectly working.

Perhaps I'm missing something...

UPDATE:

As indicated in the What's New for Geodatabases topic:

  "The recommended method to access enterprise geodatabases is to
  connect directly from ArcGIS clients. If you only use direct
  connections to your geodatabase, you do not have to install the ArcSDE
  application server."

If you need to edit settings in the ArcSDE system tables, though, you should still install the command line tools included with the application server. I don't think editing system tables directly through SQL is supported technically.


texts: QGIS: Importing geotagged photos 'Could not open' error

texts:
I am having a hard time importing geotagged photos using the 'Import geotagged photos' tool within QGIS 3.10. The photos I am importing do have GPS information within their properties, so I can't put my finger on why they are not being uploaded. This tool seems to work just fine when importing sample photos provided by my professor, but not with photos taken with my own iPhone. This message (in red) is popping up:

Can anybody help me?
In reply to HeikkiVesanto's comment - To transfer the photos from my phone to my laptop I first uploaded them to OneDrive then downloaded them from my laptop. Perhaps this is why they are in JPG. Also, the 'Import geotagged photos' tool requires photos to be in jpeg format so I don't think the image format is the issues. Also, yes my location was on. In fact, when I enter the photo's properties on my computer the longitude, latitude and altitude are all visible.
In reply to @GabrielC’s suggestion, I had tried placing the folder in my desktop and importing the photos from that location, to no avail.


texts: Replace specific raster class by surrounding classes

texts:
I have a classified landcover raster which also contains the classes "clouds" (pixelvalue 5) and "shadows" (pixelvalue 10). I want to eliminate them by replacing these two classes by the majority of the surrounding classes e.g. "forest" (value 3) or "agriculture" (value 7). I want to do this automatically for the whole raster.

It does not matter whether the clouds are replaced only by the "biggest" neighbouring class or by percentage of the "biggest" and "smaller" neighbouring classes.

Here is an example: White are clouds, black shadows, green forest and yellow agriculture.

Cloud-Masking or mosaicing was not an option in my case. So, how can I postprocess the classification and elminate these classes using QGIS? Is it possible to perform a majority filter on a specific class?


texts:
In the attribute table for a layer I have 13 Objects and if I switch on object number button. I have 16 Objects and this is wrong. I couldn´t find a way to correct this?

On the left side you see VK with 16 objects, but in the drawing there are only 13 red circles.


texts: Convert SpatialPolygonsDF boundaries to SpatialLinesDF keeping information on polygons

texts:
I have a SpatialPolygonsDataFrame and I need to convert polygon boundaries into a SpatialLinesDataFrame retaining the attributes of the 2 polygons

The structure of the SPDF looks like this

'data.frame':   11 obs. of  2 variables:
$ ID_poly  : int  1 2 3 4 5 6 7 8 9 10 ...
$ Code: chr  "33200" "32110" "33300" "32410" ...

The structure of the SLDF should be like this

'data.frame':   n obs. of  3 variables:
$ ID_line  : int  1 2 3 4 5 6 7 8 9 10 ...
$ Code1: chr  "33200" "32110" "33300" "32410" ...
$ Code2: chr  "32110" "33300" "32410" "32420" ...

where Code1 and Code2 are the values in SPDF@data$Code of the polygons that share the same boundary.

I think that

boundr <- as(SPDF, 'SpatialLinesDataFrame')

followed by an overlay (raster::intersect?) might be a solution, but I'm not sure how to properly use intersect within the same object to get the desired result (a for looping over all lines?).

This procedure has to be applied to a SPDF counting more than 10000 polygons...

This is what I'm looking for

SLDF@data
  ID_line Code1 Code2
0       1 11100 32100 # let's say, the green line
1       2 11100 41200 # same, the red line
2       3 33300 32410
3       4 32410 32420


texts: PostGIS ST_LineMerge ignoring z coordinate?

texts:
It seems that the function ST_LineMerge in PostGIS ignores the z-coordinate.

Is this a bug, or am I missing something?

Examples:

(a) This works correctly as the x, y and z-coordinates of the end points match (1, 1, 0):

SELECT ST_AsText(
        ST_LineMerge(
            ST_GeomFromText('MULTILINESTRINGZ((0 0 0,1 1 0),(1 1 0,2 2 0))')
    )
);

LINESTRING Z (0 0 0,1 1 0,2 2 0)

(b) In the example below, they last points do not match, and yet the the line segments are joined by ST_LineMerge. The result is the same as the previous example even though the z-coordinate of the first point of the second line segment is different.

SELECT ST_AsText(
        ST_LineMerge(
            ST_GeomFromText('MULTILINESTRINGZ((0 0 0,1 1 0),(1 1 9,2 2 0))')
    )
);

LINESTRING Z (0 0 0,1 1 0,2 2 0)

This is wrong surely!

Version info (on Windows 10):

SELECT version();
PostgreSQL 9.6.3, compiled by Visual C++ build 1800, 64-bit

SELECT PostGIS_full_version();
POSTGIS="2.4.3 r16312" PGSQL="96" GEOS="3.6.2-CAPI-1.10.2 4d2925d" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.2.3, released 2017/11/20" LIBXML="2.7.8" LIBJSON="0.12" LIBPROTOBUF="1.2.1" RASTER


texts: Error loading dojo file locally: "Loading failed for the  with source ..."

texts:
I'm attempting to add ArcGIS maps to a website that's on a restricted network (and therefore needs to work offline). I've been trying to get the ArcGIS ESRI maps running offline, but I'm not sure how it's supposed to be implemented. I added the arcgis_api_js folder to the root folder of my Visual Studio project (following the "AMD files hosted locally" instructions here: https://developers.arcgis.com/javascript/latest/install-and-set-up/), but when I try to reference the file mentioned, I get the error Loading failed for the <script> with source "http://localhost:57698/arcgis_js_api/library/4.18/dojo".
In my init.js and dojo.js files (referenced on the web page using the map) I set the baseURL to the address above with the intention of changing it to the IIS address used by our restricted network when live. IS there something I missed or is it not possible to run the maps locally without hosting the api library on a server at the same time?
For reference my .aspx page has this reference to the dojo file:
<script src="../arcgis_js_api/library/4.18/dojo/dojo.js">
</script>

And the dojo.js and init.js have this code replacing the boilerplate baseURL:
{async:1,baseUrl:"http://localhost:57698/arcgis_js_api/library/4.18/dojo",hasCache:{"config-selectorEngine":"lite","config-tlmSiblingOfDojo":1,"dojo-built":1,"dojo-loader":1,dom:1,"host-browser":1}

If anyone has gotten the ArcGIS maps working on a restricted network I would really appreciate your help.


texts:
Hey I'm having problems installing the dataplotly plugin from the official repository.

Konnte Erweiterung 'DataPlotly' nicht laden

ModuleNotFoundError: No module named 'plotly'
Traceback (most recent call last):
  File "/app/share/qgis/python/qgis/utils.py", line 312, in loadPlugin
    __import__(packageName)
  File "/app/share/qgis/python/qgis/utils.py", line 744, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/basti/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/DataPlotly/__init__.py", line 25, in
    from .data_plotly import DataPlotly
  File "/app/share/qgis/python/qgis/utils.py", line 744, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/basti/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/DataPlotly/data_plotly.py", line 31, in
    from DataPlotly.gui.dock import DataPlotlyDock
  File "/app/share/qgis/python/qgis/utils.py", line 744, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/basti/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/DataPlotly/gui/dock.py", line 14, in
    from DataPlotly.gui.plot_settings_widget import DataPlotlyPanelWidget
  File "/app/share/qgis/python/qgis/utils.py", line 744, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/basti/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/DataPlotly/gui/plot_settings_widget.py", line 74, in
    from DataPlotly.core.plot_factory import PlotFactory
  File "/app/share/qgis/python/qgis/utils.py", line 744, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/basti/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/DataPlotly/core/plot_factory.py", line 14, in
    import plotly
  File "/app/share/qgis/python/qgis/utils.py", line 744, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ModuleNotFoundError: No module named 'plotly'

Python-Version: 3.7.3 (default, Nov 10 2011, 15:00:00) [GCC 8.3.0]
QGIS-Version: 3.10.3-A Coruña A Coruña, 0e1f8464

Python-Pfad:
/app/share/qgis/python
/home/basto/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python
/plugins
/app/share/qgis/python/plugins
/usr/lib/python37.zip
/usr/lib/python3.7
/usr/lib/python3.7/lib-dynload
/usr/lib/python3.7/site-packages
/app/lib/python3.7/site-packages
/app/lib/python3.7/site-packages/numpy-1.14.3-py3.7-linux-x86_64.egg
/app/lib/python3.7/site-packages/MarkupSafe-1.0-py3.7-linux-x86_64.egg

I installed plotly without errors and with these commands:

pip install plotly
pip3 install plotly

Didn't had the problem on my Manjaro-Gnome machine. Guess it depends on some python path (environment variable), but can't find anything.

Can I set the path where QGIS is searching for the python-libs?


texts: Multiply or Divide certain numbers on a projection file(?) of Japan

texts:
the title seems weird though, i am dealing with Japanese shapefile and it has its own "bizarre" projection.

I have two hints to solve this issue.

it's based on espg 4301 (http://pacificprojections.spc.int/4301#)
but the coordinates is in integer(?), it has to be transformed.
for instance, the coordinates divided by (3600*1000)
e.g.) 502888725, 128444858 => 139.6913125, 35.6791272

Within QGIS, it looks like i need to make custom CRS and fix proj.4 of ESPG 4301, which looks like below (i added +to_meter = 1/3600000)

+proj=longlat +ellps=bessel +towgs84=-146.414,507.337,680.507,0,0,0,0 +no_defs +to_meter=1/3600000

but obviously, it doesn't work.

Can anyone help for this?

when i define the projection to ESPG 4301 and open the file in ArcGIS, it shows the error message of extent issue.


texts:
When I enter point id 1790 result join is only one input.
I won from list all id_gm to be shown in atribute table of layer tacka6.

Can this be done in Qgis?

Solution

https://www.youtube.com/watch?v=OowgI1b2ihw&feature=youtu.be


texts:
I'm novice in GeoServer/JavaScript.

My config:

Win 10 Pro x64 - Apache 2.14 - Tomcat 9 + deployed Geoserver 2.10.0

Apache Tomcat proxy

Tomcat 9 --- CORS

Geoserver --- Installed Vector Tiles plug-in ( geoserver-2.10-SNAPSHOT-vectortiles-plugin ) and WPS plug-in.

Trying to figure it out how to use Geoserver vector tiles so follow the instructions as written in Geoserver tutorial - http://docs.geoserver.org/latest/en/user/extensions/vectortiles/tutorial.html#publish-vector-tiles-in-geowebcache

Unfortunately, the result is a blank page with just Title "Mapbox Protobuf - vector tiles" and +/- zoom buttons ?!

What is the cause of that?

Perhaps it isn't working because Geoserver is deployed in Tomcat.

Update:

Screenshot of GEOTOOLS_DEVELOPER_LOGGING log - http://i68.tinypic.com/waiqzm.jpg

Noticed in log that BaseURL = localhost:8080/geoserver/ ,but my Geoserver is located at localhost:82/geoserver/ (!)

Maybe that's the cause??

UPDATE #2

-- regenerate tiles just in EPSG:4326 --

GEOTOOLS_DEVELOPER_LOGGING log

Request: getMap
Time = []
Filter = null
Format = application/x-protobuf;type=mapbox-vector
Buffer = 0
SRS = EPSG:4326
ScaleMethod = null
Exceptions = SE_XML
StyleBody = null
StyleVersion = null
SldVersion = null
Filters = null
Styles = [StyleImpl[ name=polygon]]
Layers = [org.geoserver.wms.MapLayerInfo@ef7f0020]
MaxFeatures = null
BgColor = java.awt.Color[r=255,g=255,b=255]
FeatureVersion = null
Transparent = true
StyleFormat = sld
Sld = null
StyleUrl = null
SldBody = null
TilesOrigin = null
ValidateSchema = false
Interpolations = []
Height = 256
Width = 256
Tiled = false
Palette = null
Bbox = SRSEnvelope[-180.0 : -90.0, 0.0 : 90.0]
RemoteOwsURL = null
Crs = GEOGCS["WGS 84",

DATUM["World Geodetic System 1984",

SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],

AUTHORITY["EPSG","6326"]],

PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],

UNIT["degree", 0.017453292519943295],

AXIS["Geodetic longitude", EAST],

AXIS["Geodetic latitude", NORTH],

AUTHORITY["EPSG","4326"]]

FormatOptions = {}

FeatureId = null

ViewParams = null

Elevation = []

StartIndex = null

Env = {}

Angle = 0.0

RemoteOwsType = null

CQLFilter = null

Version = 1.1.1

BaseUrl = http://localhost:8080/geoserver/

Request = GetMap

Get = true

RequestCharset = UTF-8

RawKvp = {FORMAT=application/x-protobuf;type=mapbox-vector, STYLES=polygon, WIDTH=256, HEIGHT=256, LAYERS=opengeo:countries, EXCEPTIONS=SE_XML, GWC_SEED_INTERCEPT=true, REQUEST=GetMap, SRS=EPSG:4326, BBOX=-180.0,0.0,-90.0,90.0, VERSION=1.1.1, SERVICE=WMS, TRANSPARENT=true}

2017-01-13 23:43:00,738 DEBUG [org.geotools.data.shapefile] - Read lock: file:/D:/GEOSERVER_DATA/data/ne_110m_admin_0_countries/ne_110m_admin_0_countries.shp by org.geotools.data.shapefile.shp.ShapefileReader

2017-01-13 23:43:00,738 TRACE [org.geotools.data.shapefile] - Shapefile Reader has been closed. Number open: 0

2017-01-13 23:43:00,738 DEBUG [org.geotools.data.shapefile] - Read lock: file:/D:/GEOSERVER_DATA/data/ne_110m_admin_0_countries/ne_110m_admin_0_countries.dbf by org.geotools.data.shapefile.dbf.IndexedDbaseFileReader

2017-01-13 23:43:00,738 TRACE [org.geotools.data.shapefile] - Dbase File Reader has been closed. Number open: 0

2017-01-13 23:43:00,738 DEBUG [org.geotools.data.shapefile] - Read lock: file:/D:/GEOSERVER_DATA/data/ne_110m_admin_0_countries/ne_110m_admin_0_countries.shx by org.geotools.data.shapefile.shp.IndexFile

2017-01-13 23:43:00,738 TRACE [org.geotools.data.shapefile] - IndexFile has been closed. Number open: 0

2017-01-13 23:43:00,741 INFO [org.geoserver.wms] -


texts:
How can I link the GUI to my processing code as shown in the image below so that it uses any number selected (between 20 and 100) in the spinSize value space?

# Perform task
if layer.type() == QgsMapLayer.RasterLayer and self.dlg.checkPixel.isChecked():

index = self.dlg.spinSize.currentIndex()
num_cell = self.dlg.spinSize.value()

processing.run('saga:resampling', {"INPUT" : output1, "KEEP_TYPE" : True, "SCALE_UP" : 0, "SCALE_DOWN" : 0, "TARGET_USER_SIZE" : num_cell, "TARGET_USER_FITS" : 1, "OUTPUT" : output2})


texts:
I want to make custom router which while calculating route costs, would take into account additional features.

Here is the use case:
Lets imagine paid motorway. Lets suppose there are 2 toll points. When router calculates route, it should, depending on user settings, find another way avoiding toll points, or inform user that his route crosses toll points here and here..

My first idea was to extend EdgeRouter2 class and override calcEdgeCost function, where I could check for given edge id in database and see if it contains toll points and act accordingly, but I want to be able to update my map and it looks like binding toll points table with edges cannot be done reliably neither by segment id nor osm_id (Is OSM_ID also part of osm2po graph?). Now I think it would be better to store toll points using lat,lon coordinates, calculate where those points are road network, and split edges at the given points.

Does this new approach make sense? What osm2po's classes should I use to accomplish this? I have a gut feeling that in this case using DefaultRouter2 instead of EdgeRouter2 would be better as it looks like it's more vertex centric - it doesn't have calcEdgeCost but it has (private) calcVertexCost function..


texts: Displaying two layers simultaneously by clicking button in OpenLayers 3?

texts:
How to display two layers simultaneously by clicking a button in openlayers3.

I want to display two layers simultaneously by clicking respective button, output should be displayed in a same map.


texts: Selecting single country and getting statistics in QGIS?

texts:
I am a new QGIS user and have the following problem:
I have a Raster-source (Map of lightpollution worldwide) and want to select single countrys and see the specific statistics concerning this country.

Is this possible and if it is, how?


texts: Using R scripts with sextante processing framework?

texts:
I want to use R scripts from the module Sextante treatment QGIS , but I have many difficulties to configure the module in Sextante . Following the procedure described in the official documentation of QGIS. I get the following error message :

"Missing dependency.This algorithm can not be run :-(
This algorithm Requires R to be run.Unfortunately , It Seems That R is not installed in your system , or it is not Correctly configured to be used from QGIS
Click here to know more about how to install and configure R to be used with QGIS "

I have completed all prerequisites needed to run the scripts:

R is installed on my machine
Sextante module and third party applications have been installed and configured. All other third-party applications work well : Saga , OTB , Grass , except R . See attachment.

I suspect an error in the path to the binary directory and scripts. Otherwise I do not see any track.

I use QGIS 2.2 from OSGeo4W stacked version on Windows 7.


texts: Does GIS Professional need CAD experience?

texts:
As a GIS Professional, one is often expected to have experience with CAD (i.e. AutoCAD, Microstation) to the extent that they can integrate CAD data into GIS workflows and datasets. However, I often see job postings in which employers are looking for candidates that have experience actually using CAD software such as AutoCAD. I would expect that drafting would be left to people with architectural and engineering backgrounds.

In people's experience, what are employers actually expecting from GIS professionals with regard to CAD experience using CAD-specific software?


texts: Setting ArcPy environment variables with Python multiprocessing?

texts:
I have a script which uses ArcPy + Python multiprocessing.

Are ArcPy environment variables shared between separate Python processes?

If I write to arcpy.env.extent from within each process am I asking for trouble?


texts: Using GROUP BY clause in arcpy.da.UpdateCursor

texts:
I need to specify a "group by" clause for UpdateCursor arcpy.da. From here I found out that it's possible to do using an optional pair of SQL prefix and postfix. But there are no examples and I'm stuck with the syntaxis problem.

Can you help me?

I've tried various modifications of "(None, 'GROUP BY "INPUT_FID"')", but still can't get it.

Here is the part of the code I'm using:

oDSFields = ('INPUT_FID','NEAR_FID','DISTANCE', popTotalField, popBlacksField)

kernel = 0
alpha = 0.005
#, "(None, 'GROUP BY "INPUT_FID"')"
with arcpy.da.UpdateCursor("outDistanceSorted", oDSFields, "(None, 'GROUP BY "INPUT_FID"')") as cursor:

    for row in cursor:

        kernel = float(math.exp(-0.5*math.pow(row[2]/alpha, 2)))


texts:
Have mercy ... with Stackexcachange I'm not yet so familiar, so this is initially somewhat opaque. I try to do my best to organize this thread as well as possible.
Recently I've started to get WFS and Feature Access working. WFS works for now. For Feature Access I need a database...

This is the "old" thread: WFS & Feature Access -> Database connection / Database support (ArcGIS 10.2) .

To make it more easy to figure out, what's going wrong, I've started from scratch.  Until now the status quo is::

I've downloaded the file "PostgreSQL_DBMS_for_Windows_922_136133.exe" from the esri Customer Portal. This package include the PostgreSQL 9.2.2 executable installation.
So...no I have a new installation of PostgreSQL 9.2.2 64Bit without PostGIS and without modifying something except of the pg_hba.conf file to enable remote access (Remote Access is working)
I've also downloaded from the esri Customer Portal the PostgreSQL_Client_Libraries for_Windows_922_136131.exe" file, which includes the Client files for database support. I've copied the 64Bit files to the \bin directory of my ArcGIS Sever installation and I've copied the 32bit files to the \bin directory of my ArcGIS Desktop installation, which is running in my personal computer.
As mentioned...the only thing, which I've modified for my psotgreSQL Database is the pg_hba.conf file to grant access via remote connections. This is working! I can access the postgreSQL database running on my server via my personal computer using pg Admin. Here is a small overview of the database:  "postgres" is the only existing database. There I haven't modified anything.

Now I want to add the Database connection to ArcCatalog. So I open ArcCatalog , go to "Add Database Connection" and there I fill in the necessary information.
The first initial problem is: I can't add the connection to the database, because I haven't got the chance to select one database. The only existing database "postgres" (look screenshot of pg Admin; the database exists!) is not selectable via the "Add Database connection" assistant . It also doesn't work to type the "postgres" manually in the "Database field". There I get the following error: Connections to system databases are not allowed. Please provide a different database for the connection."

So....the next step would be to grant access to the postgreSQL database, so I can access it via the "Add Database Connection" wizard or, maybe better?, setting up a new database. But here I don't know which settings the database needs.

I think, it would be the best to figure out step by step what I have to do. So it's easier for you and me.
Thanks a lot and I'm looking forward for some feedback :),
Dominik


texts: Iterating rasters in ModelBuilder - wildcard does not seem to work?

texts:
I've got some problem with the model below.

I want to extract rasters from sde geodatabase, use one parameter from Copy Raster, and then export it to specified folder with preserved names as an output. I need to use all operation, like on the printscreen, because this model do 'almost' what i want (maybe Iterate can be replaced by some script?).

My problem is that I want to limit my inputs - in sde.gdb I have 35000 rasters, and each of it has almost 2GB size (I know how much I can do at once after 24h, and for the whole set it will take ~ 6 days, but it would be nice to be able to leave my computer to do this only after working hours).

I want to be able to choose range (I've got unique values like ObjectID and Name which could be used to obtain it) but I don't know what I do wrong.

In Iterate Rasters I use Wildcard parameter but seems it doesn't work. Already know that for sde I need to use "%" and "_" signs instead of "*" and "?". In wildcard space I just simply wrote " OBJECTID Like '12___' " (for e.g. to be able to choose only raster with IDs from 12000 to 12999)

Anyone can me explain what I do wrong? Or maybe I need to set something more in model (like preconditions).

Maybe I do mistake and there should be some sql expression within wildcard or I need use some more tools?

Maybe I've got something wrong with my data so I uploaded small sample (4 small rasters and saved model ~5MB) my sample data. Idea will be the same for this sample and whole .sde (except changing specific signs for .sde instead of personal geodatabase I suppose)

No matter what I do, it seems like the whole process skip wildacard function and start everytime from the first ID.

Don't know if it's matter but I use ArcGIS 10.0


texts: Accessing a secure ArcGIS Server Map service from QGIS

texts:
I am trying to access a map service my organization publishes for our customers from QGIS. The service is hosted in ArcGIS for Server 10.5 and is secured.

QGIS allows me to set up a connection, but when I try to connect, I get "The query returned no layers."  The same service, if i change it to not be secured, will connect and the layers are able to be added to the map.

Has anyone else tried accessing a secured ArcGIS Server map service from QGIS? Most questions I see seem to be connecting to open services.


texts: Only read specific attribute columns of a shapefile with Geopandas / Fiona

texts:
Having geopandas installed in my Python environment, I can read a shapefile into a geodataframe with

In:
import geopandas as gpd
myShapefile = gpd.read_file(path_to_my_shapefile)
print(myShapefile)

Out:
myShapefile as a geodataframe

Unfortunately, I have some shapefiles which contain lots of attribute columns which I don't need in the end, slowing down the reading process a lot. Is there any possibility to limit the reading of the shapefile to specific attribute columns?

In regular pandas, I could use the usecols argument to the read_csv and read_table functions to limit the reading to the specified columns, e.g.

import pandas as pd
pd.read_csv(path_to_my_csv_file, usecols=['onlyThisColumn', 'andThatColumnAsWell', 'butNoOther'])

However, using usecols in geopandas' read_file method gives an error, probably because geopandas uses Fiona to read shapefiles which does not accept the argument.

  File "C:\Python34-64bit\lib\site-packages\geopandas\io\file.py", line 13, in read_file
    with fiona.open(filename, **kwargs) as f:
TypeError: open() got an unexpected keyword argument 'usecols'

Is there any other argument or way to achieve this with geopandas/Fiona?


texts: Did I really join the table successfully in QGIS?

texts:
Following a tutorial, I learn to join the vector layer "tl_2013_06_tract" with text layer "ca_tracts_pop" which contains population information. Because the identifier variable in ca_tracts_pop is text, I create a csvt file simultaneously to interpret the format. After that, the population attribute seems to be attached to vector layer successful (it appear in the new attribute table).

However, when I try to style the vector in the latter step, I can only find two original variables of  tl_2013_06_tract in the pull-down menu  in the Graduated tab.

And when opening the database manager, I could not find the new population attribute. So I wonder if I join the table successfully at the first place


texts: Adding GPS tracking data from Excel into QGIS

texts:
I have added GPS points in QGIS however all of them are the same colour. I want the points to be coloured differently for each individual.

How do I do this?


texts: How to grey-out area outside of a Polygon?

texts:
I have a Polygon on a Map which serve as "Area of Interest", how do I "grey-out" the area outside of the Polygon? (highlighted in pink)

Here's my codes:
require(["esri/Map", "esri/views/MapView", "esri/Graphic"], function (
  Map,
  MapView,
  Graphic
) {
  var map = new Map({
    basemap: "streets-navigation-vector"
    //basemap: "hybrid"
  });
  var view = new MapView({
    container: "viewDiv",
    map: map,
    center: [-118.80500,34.02700],
    zoom: 13
  });
  var polygon = {
    type: "polygon",
    rings: [
      [-118.818984489994, 34.0137559967283],
      [-118.806796597377, 34.0215816298725],
      [-118.791432890735, 34.0163883241613],
      [-118.79596686535, 34.008564864635],
      [-118.808558110679, 34.0035027131376]
    ]
  };
  var simpleFillSymbol = {
    type: "simple-line",  // autocasts as new SimpleLineSymbol()
    color: "darkblue",
    width: "2px",
    style: "long-dash"
  };
  var polylineGraphic = new Graphic({
    geometry: polygon,
    symbol: simpleFillSymbol
  });

  view.graphics.add(polylineGraphic);
});

Update:
I end up using a custom layer and create a mask using the Polygon.
https://codepen.io/kwokhou/pen/YzGwLOb?editors=0010


texts: What are the consequences of setting the SupportedImageReturnTypes to MIME in ArcGIS Server

texts:
I am currently exploring using the ExportMap function on the rest endpoint of a Cached Map Service.( I want to export only some layer within a Mapservice).

I have found that if you change the SupportedImageReturnTypes to MIME of a mxd based mapservice, you can export sublayers of a cached mapservice. (As an aside, you can't do this, if it is a msd based mapservice; even if you change SupportedImageReturnTypes from MIME+url to MIME, the exported image has all the layers. I think it directly returns the cached image.)

The ArcGIS Server help mentions that this setting

  ...Specifies whether images will be
  returned as MIME data or written to
  disk (MIME + URL). If you choose URL,
  you must have a server directory
  specified for the configuration.

I have tested this out and saw that if the setting was set to MIME+url, then a single temp png file is created in the output director, but it changes as the user pans the map & is deleted within a few seconds.

This is different from the cache, which remains on disk.

I am wondering if it is a good idea to set the setting to MIME, or would it have some consequence on performance.

Does anyone have any experience with this setting?


texts: QGIS 2.18 Windows 10, 32 bit system - where is the GRASS plugin?

texts:
I checked these question/answers already, but they did not address my issue:

GRASS plugin missing from QGIS 2.14.3 Essen in Ubuntu 14.04
How do I get grass plugin working in QGIS for Windows?

I am on a Windows 10 sytem, 64 bit. I installed QGIS 2.18 (also 2.14 and 2.16 to try to resolve the issue). The GRASS plugin is simply NOT THERE. So, when I go into the program files of QGIS 2.18 and look in the apps folder it says GRASS 7.2.0 is installed. For QGIS 2.16 the apps folder has GRASS 7.0.4. Finally, the QGIS 2.14 has GRASS 7.0.4.

I search for GRASS in the plugin console and it's missing! What am I doing wrong?

I am not a developer, so I cannot troubleshoot at the code level...


texts:
First of all, I'm completely new to qgis2threejs.

During a field trip I've measured strike/dip for some cracks in the bedrock and mapped the positions, and now I'm trying to show the results as a disk with the strike and dip for each point.

My data looks something like this (added as a delimited layer):

Point, Strike, Dip, Point, Easting, Northing, Height

1, 220, 30 NE, 1, 392430, 7302787, 9.2
...
...

When selecting "Object Type: Disk" for the point layer I can apparently choose the column with the "Dip" and "Dip direction", but there is only 4 possibilities: Fixed value, "Easting", "Northing", "Height"

Is something wrong with my data format? How do I get to choose the "Dip"-column for the dip?


texts: Moving multiple points from different layers to same location in QGIS?

texts:
I have this scenario :

These are 2 different layers, each square has points at the corner, and the goal is to have the points of each square at the same exact location with a good accuracy (if you zoom a bit they will still be at the same location).

Right now on Qgis I have to move the point from red square to the location of the point of black square, then I have to zoom and move it again (because it won't be at the same exact location), and do this again multiples times until it's good enough. I don't want to remove any point though or merge the 2 layers, I just want the points to be at the same location with a good precision.

Is there a way to merge them easily to have them at the same location ?

(and avoid this stupid thing to move a point to the location of other point, zoom, move it again...and so on)

Or maybe is it possible to select the 2 points and move them together to the same location (like that they will be at the same location)

Here is a more detailed explanation :

As you can see, all the layers are different, they don't have the same form, they can be very big or very small. In reality it's rare there are lines for the shape of layers but it was easier to draw it this way.
Each project is different, so it won't be the same position for the different layers, the form will be a bit different and so on, there won't be the same amount of layers, but it is always the same goal.

The goal is to have all the points of the layers, next to each other, be at the exact same position, except the ones which are "outside"(not close to another point)

I thought about it and I guess the easiest algorithm would be :

Go have the whole lists of points of the project
Check all the points which aren't from the same layer nearby (random number example, all the points < 50 pixels), and snap to the closest one.
If no snap, add one to counter
Go to the next point, and redo the same

At the end you will know how many points didn't snap so you can check manually why (it should be the number points at the "outside" of the map, if not, there was an issue somewhere). Do not worry if you can't do that because the snap solution is already amazing. I'll look into it later when I have time if no solution is found yet, I'll try to find one and post it here if I success in doing it but I have doubt since I never did PyQGIS code before and maybe you can't even do this idea (because for example you can't get all the points and their layer identity, or for other reasons).


texts: Renaming shapefiles with prefix while running arcpy.FeatureClassToGeodatabase?

texts:
I'm fairly new to the arcpy world. I am trying to create a script in 10.0 that will load a dataset (file geodatabase) into ArcSDE and rename the feature classes with a prefix, i.e. San_Diego_Bay -> AJ_San_Diego_Bay. My script is below and is very basic for loading shapefiles into a geodatabase. I would like to load a file geodatabase into an ArcSDE Feature Dataset. Any help would be greatly appreciated.

import arcpy
from arcpy import env

# Set environment settings
env.workspace = "Z:/Data/Test_arcpy"

# Set local variables
inFeatures = ["Atlantic_Offshore_Areas.shp", "C_C_mauritania.shp", "E_Class.shp", "San_Diego_Bay.shp"]
outLocation = "Z:/Data/Test_arcpy/test_two_arcpy.gdb"

# Execute TableToGeodatabase
arcpy.FeatureClassToGeodatabase_conversion(inFeatures, outLocation)


texts: Bluetooth my GPS signal to windows 10 tablet using QGIS

texts: Could anyone tell me how they pair their mobile phone gps to a (windows 10) tablet that can be picked up by qgis live tracking? I have tried a few apps which i believe is broadcasting the gps to my tablet (checking against google maps) but qgis is not picking up the gps information.


texts:
I've got a feature layer of certain fields (which I converted to raster), and a raster layer of different vegetation types. I now want to compute the area of each vegetation type for each field. For this I used the tabulate areas tool which did not give me any errors, but I lost some ID's. Now I read that this seems to be a general problem with big datasets.

Still I would like to ask if someone has a solution/ workaround for this that will give me the same result in the end. I'm new to python programming and have not been able to find a solution myself.

I furthermore have some entries in my resulting tables that simply do not make sense when looking at the original data. I'm worried now that the values for all entries are wrong and I cannot use them at all.

This is what my code looks like:

 inZoneData = "Field_Rast"
 zoneField = "Value"
 inClassData = "vegtypes.tif"
 classField = "Value"
 outTable = "C:/data/areatable.dbf"
 processingCellSize = arcpy.env.cellSize = "vegtypes.tif"

 TabulateArea(inZoneData, zoneField, inClassData, classField, outTable, processingCellSize)

This is a small part of the table I get where you can see that the object ID and the value do not match and in fact I have fewer OIDs than I did original values.

Edit: The datasets are originally not in the same coordinate system, but I transformed them before taking this step


texts:
OpenJump is a nice GIS, very functional, with many updates and has a beautiful rendering. But it's not very popular, and you don't see it around in discussions, in osgeo comunity, etc. And gis.stackexchange doesn't even have a tag for it.
So, what's up with OpenJump? And what keep it's development?


texts: Setting which layers are selectable in ArcMap?

texts:
I cannot figure out how to toggle selectable layers in ArcMap 10.

In 9.3, I remember being able to choose from a list which layers I wanted to be able to select features from in the map.

Let's say I'm editing two different feature classes for a building - 'FireExtinguishers' and 'EvacuationRoutes'.

While editing, I don't want to accidentally select the 'Rooms' feature class.

I want to be able to turn off the selectability of 'Rooms', but the only way I know how to change this is to right-click a layer in the Table of Contents and choose  Selection -> Make This The Only Selectable Layer.

I would like to have finer control over what layers are and are not selectable.

Does anyone know how to do this?


texts: Is there a way to order streams using the Horton (opposed to Horton-Strahler) method in ArcGIS (10.2.2)? I know the Stream Order tool allows you to use Strahler and Shreve.


texts: Gdal_polygonize: How to filter pixels above a given value (elevation)?

texts:
Given a topographic GIS raster of one country:

Given an elevation threshold n = 50 (meters)

How to polygonize pixels when z >= n ?

gdal_polygonize.py is not clear on this point. It states This utility creates vector polygons for all connected regions of pixels in the raster sharing a common pixel value. So I guess it MUST be like z = something and only this. I don't understand where to state the condition.

Current workflow details:

# download:
curl -o ETOPO1.zip 'http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/georeferenced_tiff/ETOPO1_Ice_g_geotiff.zip'
# unzip:
unzip ETOPO1.zip
# crop:
gdal_translate -projwin -005.50 051.30 10.00 041.00 ETOPO1_Ice_g_geotiff.tif crop.tif

Related to: Topography: How to convert raster into vector layers?, Get a pixel value.


texts: Loading a large osm file in osmar : why does loading time scale quadratically with file size?

texts:
I'm trying to load huge .osm files into R in order to run advanced analysis. I did not expect the loading part to be complex, but it proved difficult.

The issue here is that the time it takes to parse the XML file (.osm is xml-like) and populate an object (here an osmar one) seems to be quadratically proportional to the file size.

require(osmar)
my_file = osmsource_file("filename")
osmar_object = get_osm("bla", source = my_file)

And in case you wonder how I found that it was quadratic, here is a representation of the analysis, with files between 0 and 500 Mo.

At this point, my only chance to make it work seems to be the following : split the country in 100x100 small parts, and load them separately. Obviously, I don't really like the idea.

Why would it be so long ? How can I do so, fast ?


texts: Field calculator: generate n fields reporting n-largest values

texts:
I have a polygon shapefile and I need to generate n new fields reporting the n largest values across k variables. In this particular case n = 3, so that I would like to create three new fields: "max_area_1", "max_area_2", "max_area_3", with the highest, the second and third largest, values across variables x1, x2, ..., xk.

To make things clearer, consider this simpliefied version of my data:

x1  x2  x3  x4
10  31  73  14
 7  69  55   2
91   9  87  13

What I would like to get is:

x1  x2  x3  x4  max_area_1  max_area_2  max_area_3
10  31  73  14      73          31          14
 7  69  55   2      69          55           7
91   9  87  13      91          87          13

Using the field calculator in Python console of QGIS, I can compute the highest value using the max function. My code is:

crops = [x1, ..., xk]
processing.run("qgis:fieldcalculator", {
    "INPUT": input_layer,
    "FIELD_NAME": "max_area_1",
    "FIELD_TYPE": 0,
    "NEW_FIELD": True,
    "FORMULA": 'max(crops)',
    "OUTPUT":  output_layer})

Any idea on how I can compute, in PyQGIS, three new fields reporting the three largest values across the relevant variables?


texts: QGIS - Intersection query in a Geopackage

texts:
In QGIS 3.10.2, I wan't to calculate geometric intersections beetween two layers in a Geopackage like a view in PostgreSQL/PostGIS. The first layer is composed with features, a name attribute and a geometry.
LayerA
fid | name
----------
1   | A
2   | B
3   | C
...

The second layer is composed with 3 overlapping features, an distance attribute (originally it's a buffer) and a geometry.
LayerB
fid | distance
----------
1   | 5
2   | 20
3   | 50

For each features of LayerA I wan't to calculate intersections with LayerB by attributes (even if this instersection is null). For example :
name | distance5 | distance20 | distance50
------------------------------------------
A    | 0.005     | 0.01       | 0.05
B    | 0         | 0.015      | 0.02
C    | 0         | 0          | 0.3
...

I tryed a query with a CASE but I have intersections results only for one distance (the largest, because of overlapping I think).
SELECT
a.name,
CASE WHEN b.distance = 50 THEN ST_Area(ST_Intersection(a.geom, b.geom)) END AS distance50,
CASE WHEN b.distance = 20 THEN ST_Area(ST_Intersection(a.geom, b.geom)) END AS distance20,
CASE WHEN b.distance = 5 THEN ST_Area(ST_Intersection(a.geom, b.geom)) END AS distance5
FROM LayerA a
INNER JOIN LayerB b ON ST_Intersects(a.geom, b.geom)
GROUP BY a.name
ORDER BY a.name

An other method return me results for each feature of LayerA by LayerB attributes, but I need to pivot results to display them as column and I don't know if it's possible with a Geopackage.
SELECT
row_number() OVER () AS id,
a.name,
b.distance,
ROUND(ST_Area(ST_Intersection(a.geom, b.geom))/10000,4) as area_ha,
ST_AsText(ST_Multi(ST_Intersection(a.geom, b.geom))) as geom
FROM LayerA a, LayerB b
WHERE  a.name <> b.distance AND ST_Intersects(a.geom, b.geom)
ORDER BY a.name

What the best way to calculate area of geometric intersections between two layers, group by a layer's attribute and with overlapping features of a second layer like this ?
name | distance5 | distance20 | distance50
------------------------------------------
A    | 0.005     | 0.01       | 0.05
B    | 0         | 0.018      | 0.02
C    | 0         | 0          | 0.3
...


texts: Getting KeyError: groups from ArcGIS API for Python

texts:
I am working on ArcGIS Notebooks

from arcgis.gis import GIS
from IPython.display import display
agol_gis = GIS("https://*********.maps.arcgis.com/home/")
agolSiteAdmin = agol_gis.users.get('siteadmin')
print(agolSiteAdmin.groups)

I am getting following error

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
/opt/conda/lib/python3.6/site-packages/arcgis/gis/__init__.py in __getitem__(self, k)
   5714         try:
-> 5715             return dict.__getitem__(self, k)
   5716         except KeyError:

KeyError: 'groups'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
<ipython-input-159-b0f847f656f3> in <module>
      1 # Print all groups
----> 2 print(agolSiteAdmin.groups)

/opt/conda/lib/python3.6/site-packages/arcgis/gis/__init__.py in groups(self)
   5883     def groups(self):
   5884         """Gets a list of Group objects the current user belongs to."""
-> 5885         return [Group(self._gis, group['id']) for group in self['groups']]
   5886     #----------------------------------------------------------------------
   5887     def update_license_type(self, user_type):

/opt/conda/lib/python3.6/site-packages/arcgis/gis/__init__.py in __getitem__(self, k)
   5717             if not self._hydrated and not k.startswith('_'):
   5718                 self._hydrate()
-> 5719             return dict.__getitem__(self, k)
   5720
   5721     def __str__(self):

KeyError: 'groups'

I am using the following version

# packages in environment at /opt/conda:
#
# Name                    Version                   Build  Channel
ipython                   7.5.0            py36h39e3cac_0    defaults
ipython_genutils          0.2.0                    py36_0    defaults


texts: How to convert a shapefile into a graph in which use Dijkstra

texts:
I am working in a simulation using NetLogo, the idea is to make the agents interact inside the roads of a simplified city, in order to  do that I need to use some algorithms to find the fastest path between two points in the city (for example Dijkstra or A*), and to do that I need to work with a graph of the city (with nodes and vertices) but I don't know how to convert the shapefile into a graph.

My idea was to get all the intersections between two streets and start building the graph connecting those nodes. I have tried networkx for python but I had some problems installing the library gdal. Currently I'm using the s2g library (here) to do it, but for example for a shp file like this

I always get this:

INFO:root:Validating pair-wise line connections of raw shapefiles (total 3 lines)
100% (3 of 3) |##########################################################################################################| Elapsed Time: 0:00:00 Time: 0:00:00
Major components statistics:
    Total components: 0
    Component size: max 0, median 0.0, min 0, average 0.0
    Top comp. sizes: 0

I suppose that is an empty graph so I cannot work with it. Perhaps there is a better approach to my problem or an easy tool to convert shp to graphs. I mac user so I have no access to spatialite.


texts: How to adjust rows and columns, retaining cell sizes in arcmap?

texts:
The Situation: I have a vector data map with all the countries in the World, including a unique landcode for every country. I need to convert it into a raster map, so I can eventually analyse (with Matlab) drought indicators for that country.

The Problem: The precipitation data I use has 720 columns and 360 rows. I want the raster map of world countries to be in the exact same amount of columns and rows. BUT the raster map also needs to have a cellsize of 0.5! Right now, when I use the vector to raster conversion tool, if I use cellsize = 0.5, the columns and rows become 720x287.

Question: How can I create a raster map of 720 columns and 360 rows, while setting cellsize at 0.5?

Additional info: Using ArcMap 10.4.1
Using Matlab R2015b
Raster is .TIFF
Vector is .shp


texts: How do I redefine the projection of a DEM?

texts: I have a DEM that was exported from AutoCAD civil 3D 2013 but it's in the wrong projection. I'm using an area specific projection that was not provided as an option in CAD DEM export dialog box. I choose an other projection to finish the export and now I want to associate the proper projection to the DEM. Is there a way to redefine the projection of the DEM in ArcGIS 10.1 to the proper one?


texts: How to place labels at regular intervals orthogonally to a line feature?

texts:
i want to display the current chainage of a line in a certain interval (say every 100m). I know this can be done for instance by the plungin Qchainage or LRS. But the labels have to be displayed orthogonally to the current line bearing, which none of these tools do provide fields for. Considering the advanced labelling options i know i can set labels (even with flexible bearing option) at certain intervals but I did not find an option (expression) to display the respective position of a line.

Is there a straight forward way to solve an issue like I delineated?


texts:
I am new to SQL. I am trying to execute some queries on a database using PostgreSQL with the PostGIS extension. I cannot understand why the execution of the following query takes about 5 mins

select st_buffer(l.geom, 0.015), st_buffer(s.geom, 0.015)
from   "road" as l,  "building" as s

while the execution of these queries

1) select st_buffer(l.geom, 0.015)
   from "road" as l

2) select st_buffer(s.geom, 0.015)
    from  "building" as s

take about 300 ms for each query.

I have noticed that the execution of the query below takes a long time too.

select st_buffer(l.geom, 0.015)
from   "road" as l,  "building" as s

Moreover, I want to execute the ST_Intersection between the two buffer resulting from the queries. Something like this:

select st_intersection(st_buffer(l.geom, 0.015), st_buffer(s.geom, 0.015))
from   "road" as l,  "building" as s

What is the correct way to do that?


texts: Annotation feature class interfering with Maplex labels?

texts:
Because an annotation feature class is static and separate from the dynamics of a labeling engine, I'm having trouble understanding why my street annotation is altering my labeled features.

More specifically, I have milepost point features that use an expression to essentially place the numbers vertically, like you see in real space when driving down the highway. And it works just fine when the street annotation layer is turned off. But when I turn it back on, some of the labels ignore the expression and revert back to a standard horizontal, left-to-right labeling of the numbers (see related images and note milepost 298 at bottom of map page).

How do I fix this so that the mileposts always (i.e. hard constraint) label according to the expression?


texts: Explaining what GIS is to 11 year old kid?

texts:
Well, since it was one of the highest rated proposed questions, hasn't been asked yet, and I would love to know the answer, I thought I'd ask it.

How do I explain what GIS is for the 11 year old kid?


texts: mosaic dataset and overview failed in arcgis

texts:
I m working with arcgis 10.
i m trying to optimize the storage and the visualization of heavy raster.
After importing raster in file geodatabase i put them in a mosaic dataset and i generate overview but:

Some raster are missing, and some are cut.
Overview are shifted

However i m sure that my input raster are in good projection, with stastistics calculated and pyramids.
Any help?


texts: What should I pass to ILayerFactory.Create for a WMSLayerFactory?

texts: I was expecting ILayerFactory.CanCreate to return true for an IWMSLayerDescription.  I recursively traverse all the descriptions in a WMS connection but CanCreate always returns false.


texts: Clustering trajectories (GPS data of (x,y) points) and mining the data

texts:
I've got 2 questions on analyzing a GPS dataset.

1) Extracting trajectories I have a huge database of recorded GPS coordinates of the form (latitude, longitude, date-time). According to date-time values of consecutive records, I'm trying to extract all trajectories/paths followed by the person. For instance; say from time M, the (x,y) pairs are continuously changing up until time N. After N, the change in (x,y) pairs decrease, at which point I conclude that the path taken from time M to N can be called a trajectory. Is that a decent approach to follow when extracting trajectories? Are there any well-known approaches/methods/algorithms you can suggest? Are there any data structures or formats you would like to suggest me to maintain those points in an efficient manner? Perhaps, for each trajectory, figuring out the velocity and acceleration would be useful?

2) Mining the trajectories Once I have all the trajectories followed/paths taken, how can I compare/cluster them? I would like to know if the start or end points are similar, then how do the intermediate paths compare?

How do I compare the 2 paths/routes and conclude if they are similar or not. Furthermore; how do I cluster similar paths together?

I would highly appreciate it if you can point me to a research or something similar on this matter.

The development will be in Python, but all kinds of library suggestions are welcome.

I'm opening the exact same question https://stackoverflow.com/questions/4910510/comparing-clustering-trajectories-gps-data-of-x-y-points-and-mining-the-data in StackOverflow. Thought I'd get more answers here...


texts: Displaying map from ArcGIS REST API using Python

texts: Is there any way to display an ArcGIS map from a standalone ArcGIS Server URL ie (https://xxx.xx.xxx/acgis/rest/services/MapServices/xxxxxmap/MapServer/) using Python? ideally ipyleaflet.


texts: How to reference vector based data to raster data that already has spatial reference in MAPublisher?

texts: I have been able to load in raster data into MAPublisher and properly registered it, however I have been having trouble with a further step. I am trying to reference vector data that has no GCS or projection to the  the raster data's spatial information. The vector data is in pdf format, and I would like to be able to reference it and export it as an ESRI supported format. I have been able to size the vector data to match the raster data, but every time I try to register it to the raster data, it makes it incredibly large and goes off the canvas. Is what I want to do possible?


texts: Converting OSM.PBF to Geopackage in CLI?

texts:
I would like to refresh this tutorial for my purposes - background vector map in one file - just to load in QGIS.

My idea is a script doing this:

Import directly downloaded Openstreetmap PBF file (country PBF from Geofabrik) or cropped (selected region etc.) using ogr2ogr directly to GeoPackage or some temporary working file.
simplify the attribute data - if I need just a basic background map with labels I do not need many attribute data fields and can merge others.

Should I use some temporary file (like Sqlite) for the attribute data manipulation or can I do this with Geopackage directly using set of CLI commands?


texts: Updating folium changed the outlook of Popup box

texts:
Recently I updated folium and noticed some changes to the popup box. With the new update the popup box seem to have shrinked and the text is coming in separate lines, which should have appear in the same line as before. There was no changes been made with the code. How can I change the popup box look like the previous one, i.e., appear in the box without breaking the text?

Code for the Marker Popup:
fgc.add_child(folium.Marker(location=[lt, ln], popup= "
 Thana :&nbsp" + di +"
"+ "
Cases Total: &nbsp: "+str(ca)+ " person "+ "
" + "Cases 24 hours : "+ str(da)+ " person "+"
"+"Cases 7 days: "+str(we)+ " person "+"
Neighbouhood affected :  "+str(ne),


texts: ESRI geoprocessing service with Javascript API 3.5

texts:
I am having issues getting a geoprocessing service to run.  The service is a simple Python script that converts a table to CSV file.  It runs fine in ArcMap. I can publish the service with the input as a constant and get it to run from Javascript. The problem is when I try to run the service with Feature Set created from a Feature Layer in Javascript. I can not get this to work!

My geoprocessing service is here.

function exportAll(){
    var query = new esri.tasks.Query();
    query.where = "1 = 1";
    query.outFields = ["*"];
    //parcelsFS is my Feature Service Layer
    parcelsFS.queryFeatures(query, function (features) {
    var inputFeatures = new esri.tasks.FeatureSet();
    inputFeatures = features.features;
    var params= { "Table_Name":inputFeatures};
    gpTask.submitJob(params, completeCallback , statusCallback,function(error){
      //alert(error);
      //esri.hide(loading);
        });
    });
}

function completeCallback(jobInfo){
    if(jobInfo.jobStatus !== "esriJobFailed"){
      gpTask.getResultData(jobInfo.jobId,"Output_CSV", downloadFile);
    }
  }
  function statusCallback(jobInfo) {
    var status = jobInfo.jobStatus;
    if(status === "esriJobFailed"){
      alert(status);
      //esri.hide(loading);
    }
    else if (status === "esriJobSucceeded"){
      //esri.hide(loading);
    }
  }
  function downloadFile(outputFile){
    map.graphics.clear();
    var theurl = outputFile.value.url;
        console.log(theurl);
    window.location = theurl;
  }


texts:
I want to show tooltip for my dojo button, I have tried the following code but getting error: Uncaught SyntaxError: Unexpected token (. Please solve this issue.

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no"/>
    <title>Simple Map</title>
    <link rel="stylesheet" href="http://js.arcgis.com/3.14/esri/css/esri.css">

 <button id="Info" data-dojo-type="dijit/form/Button">Draw PolyLine</button>


texts: The proportion of population covered by created isochrones in QGIS

texts: I created isochrones using Iso-Area as polygons algorithm is QGIS. I would like to calculate the proportion of the population covered by the created isochrones. Is there any tool in QGIS to do this?


texts: How can I create contours from GeoTIFF and Python (GDAL, Rasterio, etc?) into shapefile with prj file

texts:
I'm using a few reference StackExchange/Overflow posts here - and manual, but not coming up with what I'm after.

I have a GeoTIFF that I'm importing, and I need to develop out contours at regular intervals.

I've repurposed code from this post, but the two issues I'm having are:
- I can't seem to specify regular intervals, rather I have to load my intervals manually into an array
- The shapefile is created without a prj file, meaning it's not spatially referenced.

My code so far:

from osgeo import gdal
from osgeo.gdalconst import *
from numpy import *
from osgeo import ogr

#make contours from the raster
print('Making Contours')
indataset1 = gdal.Open(os.sep.join([output_dir, os.path.basename(sabre_out_dir)]) + 'mdROS.tif', GA_ReadOnly)
in1 = indataset1.GetRasterBand(1)
dst_filename = os.sep.join([output_dir, os.path.basename(sabre_out_dir)]) + 'mdROS_contour.shp'
ogr_ds = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(dst_filename)
contour_shp = ogr_ds.CreateLayer('contour')

field_defn = ogr.FieldDefn("ID", ogr.OFTInteger)
contour_shp.CreateField(field_defn)
field_defn = ogr.FieldDefn("mdROS", ogr.OFTReal)
contour_shp.CreateField(field_defn)

#Generate Contourlines - Instead of manually specifying all intervals I'd like to say "every X interval"
gdal.ContourGenerate(in1, 0, 0, [1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000,9000,10000], 0, 0, contour_shp, 0, 1)
ogr_ds.Destroy()


texts: Workflow automation QGIS - populating attribute table

texts:
Is there a tool in QGIS which populates fields within the attribute table?
I know I can edit it manually but what I want is to include that tool in the workflow automation process.
I want to label the output rather than putting in on the key of the map. In order to label I need to have a field where I can add the label. Since this is part of a larger process in the Graphical Modeler I need to automate it rather than adding the field and name manually. What I'm interested in is if there is a tool or Expression in the Field Calculator that can return the name of the input layer and add it to the attribute table of the new/ output layer of the Field Calculator tool


texts: GRASS GIS: Getting know about stream power index (SPI) computation by r.watershed function

texts:
I need to know exact way how stream power index (SPI) is computed by module r.watershed. On the help page is that SPI is α * tan(β), where:

α is cumulative upslope area draining through a point per unit contour length

and

tan(β) is the local slope angle.

So I tried to make an example digital elevation model (DEM) to confirm my computation "by hand" and numbers computed by r.watershed. So I made 1 000 meters long and 200 meters width DEM, where upper pixels has 150 meters above sea level and lower pixel has 100 m.a.s.l. So slope at this example dem is 5% or 2.862°. Pixel resolution is 5 x 5 meters.

As α I consider flow acummulation (also computed by r.watershed) and I want to compute SPI value for very last pixel in flow direction (pixel with 100 m.a.s.l.). Flow accumulation number was 198,5. So, by pixel resolution 25 square meters, cumulative upslope area is 198,5 * 25 = 4962.5 square meters.

As β I take slope in degrees, so 2,862°. Therefore, tan(β) is 0,0499.

By the equation on the help page, SPI = 4962.5 * 0,0499 = 148.875.

But value from SPI raster is 16.51.

I tried to find some explanation, but I can't find anything. I also try to find something in source code, but its pretty hard for me.

Do somebody know where can I find some additional info about how is SPI computed? I assumed that point of my question is hidden in α value, but I can't figured it out.


texts: How to set drawing style for QgsRasterLayer in QGIS 2?

texts:
I am using PyQGIS.
In QGIS1, I can use:

rlayer.setDrawingStyle ( QgsRasterLayer.PalettedColor )

but In QGIS2, I got some errors.

I follow pyqgis developer cookbook, put some demo scripts in qgis python console without success.

rlayer.drawingStyle() #got AttributeError
rlayer.setDrawingStyle(QgsRasterLayer.SingleBandPseudoColor)
#got AttributeError from QgsRasterLayer, change it to QgsRaster
rlayer.setDrawingStyle(QgsRaster.SingleBandPseudoColor) #got TypeError


texts:
I am trying to join a shapefile and a CSV file. I have tried Arcpy.AddJoin_management... the error message I receive is that there is no field name for the CSV file.

I tried to index it and same error message shows up.

import arcpy
from arcpy import env

env.workspace = "E:\county"
in_data = "county.shp"
out_data = "E:\county\County2.shp"
data_type = ""

arcpy.Copy_management(in_data, out_data, data_type)
arcpy.env.overwriteOutput = True
arcpy.MakeFeatureLayer_management ("E:\\county\\County2.shp")
arcpy.MakeTableView_management ("E:\county\data.csv")

LayerName = "County_Layer"
Field = "NAME"
TableName = "data_View"
Field2 = "County"

arcpy.AddJoin_management (LayerName, Field, TableName, Field2)

Traceback (most recent call last):
  File "<pyshell#15>", line 1, in <module>
    arcpy.AddJoin_management (LayerName, Field, TableName, Field2)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\management.py", line 6066, in AddJoin
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000728: Field County does not exist within table
Failed to execute (AddJoin).


texts: FME: Can I add multiple feature types to the SQL Server Spatial reader?

texts:
The title's it, really.

With .TAB file reader, you can add multiple layers to the reader and use the 'merge feature type' check box with an '*' to read all the layers listed in the source files, which can then be separated out again on featuretype name.

With the SQL Server Spatial reader, you can specify many layers in the database selection, but it seems that the 'merge feature type' option doesn't work.

Is this a bug, or expected behavior, or am I missing something?

EDIT: I have found that I can do this sometimes... so the above may be incorrect. I'll need to test some more.


texts: How to fix symbols of layers when I copy .mxd to another pc?

texts: I've created a personal geodatabase with feature classes and I specified symbols for them in an mxd project. I activated "store relative pathnames to data sources" in the Map Document Properties. When I sent the project and database to another computer, the symbols of feature classes were changed. Is it possible to fix this?


texts: Extract areas from DEM between 100 meter contours using r.report or similar tools

texts:
I'm trying to extract the total area within the same altitude range in a DEM. As example, I'd like to know the total area between 800-900 m asl of a DEM.
Using r.report I've just got this:

How can I set up the range values instead of a 'random' ranges generated by the gis tool?

I'd like to use this in the same way to calculate slope ranges etc.


texts: Showing only selected features from attribute table on layer using QGIS?

texts:
One of the layers has around 5000 attribute entries, based on a CSV file of species records. All the features are of course shown on the map layer. How can I select a specific class of features so that only they appear on a map, so that for example a distribution map for a species can be generated. I have tried "Advanced Filter (Expression)" so that the Attribute Table shows only the required species, but it still shows all the attributes on the map layer.

I am using QGIS 2.10.1.


texts: Labeling selected items out of a column in an attribute table in QGIS

texts:
I have a layer of roads from OSM.
I have selected 166 "Motorway or highway" road type out of 3,986 rows.

I need to label only this road type without deleting any data.
When I am typing  "TYPE" in the "Label with" window I have all the road types labeled.
How can I label only this selected rows?


texts: When I add dgn file in arcmap, nothing shows up on my map. what should I do with the dgn file so I can see the feature on arcmap?


texts: Iterating a model that contains values too low

texts:
Beginner Python programmer here,

I have exported my Model Builder model from ArcGIS to Python. What I am trying to do is provide some kind of clause saying "if the data selected has less than 2 records/features, skip it and go to the next iteration."

I have found information on SearchCursors and CurrentStates, but I am look at the entire feature. Not just the one field.

Any help?


texts: Trouble getting lat/lons from point layer

texts:
I am following this post about how to get lat/lon information from a point file, but without success.  Does this metadata posted below provide information that changes the prescribed method for getting lat/lon for each point?

In layer spatial reference system units : xMin,yMin 1091130.77,1813891.89 : xMax,yMax 1205198.85,1951669.02
Layer Spatial Reference System:
+proj=tmerc +lat_0=36.66666666666666 +lon_0=-88.33333333333333 +k=0.9999749999999999 +x_0=300000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs

I have taken the following steps:

Loaded a ESRI shapefile into QGIS with one polygon layer.
Extracted nodes to create a point layer. Vector -> Geometry Tools -> Add/Extract Nodes. I clicked the option to create new shapefile, then clicked the option to add shapefile to the current project.
Right clicked on new point layer, selected 'Save As', and specified target CRS as ESPG:4326, and added GEOMETRY=AS_XY in the layer OGR Creation options. This created point layer 2.
Clicked on point layer 2, then ran Vector -> Geometry Tools -> Export/Add Geometry Columns.

I ran this entire process four times, once with no layer OGR creation options, and once with the Geometry=AS_XY in the layer OGR creation option.  I also varied step 4: once with target CRS 'layer' in step 4, once with CRS 'project' in step 4.

I also added a step 5, Right click on point layer 2, then save as CSV file.

When complete, I right click on point layer 2, click 'Open Attribute Table', and see new XY columns with numbers that are definitely not lat/lons.  Also no lat/lons in the CSV file.


texts: How to view a part of the map in natural size using ArcGIS for Desktop?

texts:
I have a map project in ArcGIS 10.3. My project page size is custom 1500 * 1000 mm.

How can I cut a part of this map to print it on A4 format paper and see this part with a natural resolution (same as it will be in 1500 * 1000 mm size)?

EDIT: Added picture for clearing what I want.


texts: Does anyone know what format of raster images works with Orux Maps? I tried to export maps acc to Orux Maps specs from ArcGIS 10.1 in every possible format (Geotiff, GRID, jpeg2000...), loaded with a .tfw file and it always displayed a colorless tile.


texts:
In my dataset there are two points that are identical in every way. I'm looking for a way to run some kind of query or functional attribute to make it such that there is only one point at this location.

Reason is that if I try and label this data set there will be two identical values and because they are 'Read-Only' types I cannot simply delete one of the labels.

I am working in GeoMedia Professional 2014


texts: A count of points within a set distance using updatecursor

texts:
I have a script that calculates the average distance between points. What I want is to rewrite it so I can input the point layer, input a distance, and the output field will tell me a count of the number of points at that distance. Can anyone help? so basically I just want a script that will tell me a count of the number of points from a point at a certain distance.

#Import standard library modules
import sys, string, arcpy

    inputFC = arcpy.GetParameterAsText(0)# Input Point Layer
    inputdistance = arcpy.GetParameterAsText(1) #input distance
    outputfield = arcpy.GetParameterAsText(1) #output field

# Read the point data into list.
#
    xs = []
    ys = []
 rowCursor = arcpy.SearchCursor(inputFC)
 for row in rowCursor:
    geom = row.getValue (properties.shapeFieldName)
    cen = geom.centroid
    xs.append(cen.X)
    ys.append(cen.Y)
 del rowCursor

 # compute distance
 sumdist = 0.0
 count = 0.0
 for i in range(0,len(xs)):
    for j in range(i,len(xs)):
        if (i <> j):
           dist = ((xs[i] -xs[j])**2 + (ys[i] - ys[j])**2)**0.5
           sumdist = sumdist + dist
           count = count + 1
  print "The average inter-point distance is " + str(sumdist/count)
  arcpy.AddMessage("The average inter-point distance is " + str(sumdist/count))

With the following script. Ive tried to run it as  Distance.add_nearby_points_count_field('gauges.shp', 100)

 import sys
 import arcpy
 import math

def add_nearby_points_count_field(inputFC, inputdistance):

# Read the point data into list.
coord_pairs = []
desc = arcpy.Describe(inputFC)
rowCursor = arcpy.SearchCursor(inputFC)
for row in rowCursor:
    geom = row.getValue(desc.shapeFieldName)
    coord_pairs.append(geom.centroid)

# Add a field called 'cnt' and calculate nearby points
arcpy.AddField_management(inputFC, 'cnt', 'LONG') #i also tried cutting   out this part and creating the "cnt" field manually
rowCursor = arcpy.UpdateCursor(inputFC)
for row in rowCursor:
    geom = row.getValue(desc.shapeFieldName)
    from_point = geom.centroid
    near_pts = 0
    for to_point in coord_pairs:
        distance = math.sqrt(pow((to_point.X - from_point.X), 2) + pow((to_point.Y - from_point.Y), 2))
        if distance <= inputdistance:
            near_pts += 1
    row.cnt = near_pts - 1 # Subtract 1 to remove the measurement to itself
    rowCursor.updateRow(row)
del rowCursor

if __name__ == '__main__':
    inputFC ='guages.shp' #originally its a featureclass but I converted it to a shapefile to try both methods
    inputdistance = 100
    add_nearby_points_count_field(inputFC, inputdistance)


texts: How to create the fuel models & live fuel moisture rasters for Wildfire simulation?

texts: Has anybody done wildfire modelling in GRASS?  Where could I obtain / how do I create the fuel models and live fuel moisture rasters to determine rate of spread (r.ros)?