texts: Collecting values from iterated intersections in ModelBuilder?

texts:
I have a points feature class which I am intersecting with a number of polygon feature classes, one at a time (i.e. points AND polygons1, points AND polygons2, etc...).

Each time I do an intersection I note the number of points captured by the polygons layer I used. I have a modelbuilder model that goes through the first part of this process iteratively.

How can I collect the number of points intersected, relative to the polygon layer I used, as a list/table?

I'm trying to use the Collect Values tool to tabulate the number of points per intersection but I don't understand where the information is going. Should I be using a different tool?

I noticed that Get Count may be the tool I need to get the number of records geenrated after each intersection. However, I don't know how to output this information from successive iterations into different rows of a table.


texts: PostGIS: Running multiple intersections within a single SQL query?

texts:
PostgreSQL 11 + PostGIS, on Windows 10.

I'm working with a PostGIS-enabled database.
I want to UPDATE a large polygon table (100,000+ rows) WHERE the geometry of that host table intersects the geometry of another table (table a) AND another table (table b). 'Table a' also has a WHERE field = value attribute filter placed on it.

Using this statement, it takes an indefinte amount of time, seemingly never to complete (I even left it overnight and it still didn't finish):

UPDATE schema1.polygontable poly
SET ref_id = '5z'
FROM schema2.tablea a, schema3.tableb b
WHERE (ST_DWITHIN(a.geom, poly.geom,0) AND (a.type = 'Ground')) OR (ST_DWITHIN(b.geom, poly.geom,0));

However, If I split the query at the OR clause, the now separate 2 functions take 5 seconds and 1 second respectively:

UPDATE schema1.polygontable poly
SET ref_id = '5z'
FROM schema2.tablea a
WHERE ST_DWITHIN(a.geom, poly.geom,0) AND (a.type = 'Ground');

and

UPDATE schema1.polygontable poly
SET ref_id = '5z'
FROM schema3.tableb b
WHERE ST_DWITHIN(b.geom, poly.geom,0);

Nb: all three tables have a spatial index each. And they are all hosted on an Amazon AWS instance.

So, why does my query freeze/hang when I combine the functions? My syntax looks fine I think


texts: After export the clipped Landsat image, it's distorted

texts:
I clipped square(polygon) of NDVI and Aster GED dataset and select the scale of 100 and set CRS to native projection and Exported to drive. But when I am Opening this dataset in QGIS it's showing image size 263X49 and resolution in degrees (which is not equal to 100m). Why is my image not showing equal number of rows and columns (clip-square) and in some cases image is completely distorted?
Here is the link for my code. https://code.earthengine.google.com/17caeb1658b9b15c2177077d0731f02e?accept_repo=EE201


texts: Is it possible to intersect multiple polygon shapefiles to identify areas where two or more layers overlap?

texts:
I have 19 polygon shapefile layers, with some occasionally overlapping with one another. As a final output, I would like a shapefile that includes only the overlapping areas between 2 or more of these 19 layers (and if it could keep their respective attributes, that would be fantastic).

The Intersect tool provides an ideal output, yet if I add all 19 layers as the inputs, it searches for areas where all 19 intersect. The only remedy I can currently imagine is doing pairwise intersections, which would mean I'd have to do it 171 times, then merge them all into one shapefile layer.

I cannot seem to find a method for what I want to do that doesn't involve painful pairwise analyses.


texts:
I am unable to turn on the labels within my Python script:

import arcpy
mxd = arcpy.mapping.MapDocument(r"Mypathhere") #Map document reference
layer = arcpy.mapping.ListLayers(mxd, "Layername")[0] #Indexing list for 1st layer
if layer.supports("LABELCLASSES"):
    for lblclass in layer.labelClasses:
        lblclass.showClassLabels = True
arcpy.RefreshActiveView()
mxd.save()
del mxd

If I need to enable labeling in the label manager first, can that be done with ArcPy?

The script does not break. I do not get any errors. The points just do not get labeled.


texts: Coverage Turn Table to Turn Feature Class in ArcGIS Pro

texts:
I have data about streets in coverage format (arc, node and turntable). I would like to transform that data to ArcGIS Pro and use Network Analyst for calculating navigation.
Steps, that I used:

I created new gdb and added to it arc and node (I used import feature class)
I used function: Turn table to Turn Feature class and got Turn Feature Class
I created network dataset
As Source Feature Classes I used: arc, node and turn
When I want to calculate route, my network dataset does not recognize forbidden turns

How can I get the network to recognize forbidden turns?


texts: How can I calculate the amount of memory required to run osm2pgrouting?

texts:
I have read that osm2pgrouting loads all the data into memory, which makes it faster.

On the other side the size of the .osm file which it can parse is limited.

How can I exactly know how much memory would I need to parse for example hamburg.osm into Postgres with osm2pgrouting? (hamburg.osm is 17 mb zipped and more than 100 mb unzipped).

Is that possible to calculate the required memory size at least approximately?

P.S. I have tried to parse a 300 Mb .osm file on a Linux machine with 512 Mb of RAM. It got to "Classes table created" and after that has not returned any "Finished" message. It also created 3 tables in the database - and populated the Nodes table with about 1.5 million nodes. (I have followed the the recent pgRouting workshop).


texts:
I am using geoprocessing tool Intersect in Arcobjects. How to pass two input feature classes in to in_features option in INTERSECT tool could you please tell me how to pass.

Geoprocessor gp = new Geoprocessor();
ESRI.ArcGIS.AnalysisTools.Intersect Inter = new ESRI.ArcGIS.AnalysisTools.Intersect();
Inter.in_features = ????????????????????????????????
Inter.out_feature_class = "D:\\Temp\\Inter.shp";
gp.OverwriteOutput = true;
gp.Execute(Inter, null);


texts: HTML form to load points on Google Maps using API?

texts:
I am creating a web page with a Google Map embedded in it.  I want to create a form to allow a user to type which London tube zone they wish to see and, when submitted, this will map the tube stations for the selected zone.

I have the map embedded where I want it to be and centred with the zoom I need.  I can hard code the map to display one of the zones using:

downloadUrl("ldntube_xml.php?zone=1", function(data) {

But I want the user to be able to select the zone they want.

For the form so far I have:

<form action="ldntube_xml.php" method="get">
    Type Zone Number:

    <input type="text" name="zone" value="" >

    <input type="submit" value="Submit">
</form>

This obviously takes the user away from the web page to list the stations in the zone selected.

How can I make the form and the map display talk to each other?
Can I accomplish this using an onClick"myFunction()"?


texts:
I'm trying to display point data as diagrams with labes next them. I'm getting the diagrams and labels to work separately, but they won't work together. When both are turned on, the labels disappear and as soon as I turn off the diagrams, the labels reappear.

Does anyone know how to show them together?


texts: OOM uploading a large amount of .tif - aerial imagery ot PostGIS with python

texts:
I'm trying to upload our big dataset of aerial imagery on our server via python and psycopg2
The function I use for this task is:
def load_into_PostGIS(connection, image_path):
    with open(image_path, 'rb') as f:
        with connection:
            with connection.cursor() as cursor:
                cursor.execute("INSERT INTO remote_sensing.aerial_imagery(rast) VALUES (ST_FromGDALRaster(%s))", (f.read(),))

Uploading up to 3 of these around 5 MB. sized .tiff-files works with no issues. After that, it's impossible to open them in qgis and actually raises an OOM-Error after some while when I run "load_into_PostGIS" with .apply on my pandas Dataframe containing all the files paths.
What am I missing here that uploads those file obviously not the way I need them?
Edit:
Instead of using the funky python Insert statement I use raster2pgsql:
    cmd = f"""
            set PGPASSWORD={pw}&&SET PATH=%PATH%;{bin_path}&& "raster2pgsql.exe" -I -C "{main_path}\dop*.tif" -F {schema}.dop20_rgb_utm32_{folder} | psql -U {user} -d {dbname} -h {host} -p {port}
           """


texts:
I am attempting to run Extract by Mask through a loop using Search Cursor, where the Cursor is defined by the row "FID" in a table from a grid shapefile. The code is as follows:

import arcpy, os
from arcpy import env

outws = r"C:/path2/"
raster = r"C:/path1/raster.tif"
fc = r"C:/path1/shapegrid.shp"

arcpy.CheckOutExtension("Spatial")
count = 0
with arcpy.da.SearchCursor(fc, ["FID"]) as cursor:
    for row in cursor:
        out = os.path.join(outws, "demgrid" + str(count) + ".tif")
        ebm = arcpy.sa.ExtractByMask (raster, row)
        ebm.save(out)
        count = count + 1

arcpy.CheckInExtension("Spatial")

sys.exit()

The error that first pops up is on the 17th line: ebm = arcpy.sa.ExtractByMask (raster, row)

Traceback (most recent call last):
  File "grid_extract.py", line 17, in <module>
    ebm = arcpy.sa.ExtractByMask (raster, row)

My speculation is that either I've called the wrong module, or that I've not defined a specific one I need for this task.


texts: Classifying an image based on the EVI values of vegetation

texts:
I am very new at GEE, and I’m trying to classify an image based on the EVI range of various vegetation classes. For example, Bamboo may have a minimum threshold of 0.6712 and maximum of 0.776 EVI. This way I have 10 classes. I obtained the EVI band from my Landsat image.
How do I classify the pixels and assign them values from 1-10 for the 10 classes?
I tried using ee.Algorithms.If method:
    var EVIband = ee.Image(CFL598.select('EVI'));
    print("evi", EVIband);
    Map.addLayer(EVIband, {min:-0.10, max:0.8, bands:['EVI']}, 'EVI Band');

    var ten = ee.Image(ee.Algorithms.If((EVIband.gte(0.6712).and(EVIband.lte(0.776))), EVIband.mask(ee.Number(10))));
    print("bamboo", ten);
    Map.addLayer(ten);

I was planning to do this for the rest of the 9 classes and use addBands to create a whole image, but this did not work.


texts: Load a shapefile in PostgreSQL database with PHP

texts:
I am trying to load a shapefile in a PostGIS database using the exec() command. I have read the format I have to follow from here and I checked several similar questions from the forum.

Based on the above I wrote this script:

  $command = 'shp2pgsql -s 2100 -d /home/enomix/www/EsoterikoFinal/Maps/shps/nodes.shp  | psql -h localhost -p 5432 -d shapeDb -U postgres -W xxx';
  $output = exec($command);
  print_r($output);

What I want is to load the shapefile into a newly created table which is in a database called 'shapeDb'. Every time I load a shapefile I want the previous table to be dropped and a new table to be created (that's why I use the -d option). I have defined a full path for my shapefile (should I define a relative one like: "shps/nodes"?)

When I run the script I get nothing back as a response although I have a:

ini_set('display_errors', 1);

Its the first time I do something like this. Any idea what I am doing wrong?

EDITED

Running the above command on the terminal it works fine. But I still can not figure out how to drop the previous table and create a new one each time I run the query. And also why the exec command doesn't work!

EDITED II

Now I made it work. The -W xxx was causing problems. I don't know why though (its just the password). I run my script like this:

$command = 'shp2pgsql -s 2100 /home/enomix/www/EsoterikoFinal/Maps/shps/nodes.shp  | psql -h localhost -p 5432 -d esoteriko -U postgres';
exec($command,$out,$ret);

BUT still I can not figure out how to drop the table each time and create a new one. ANY IDEAS?

EDITED III

This works now:

$command = 'shp2pgsql -s 2100 -d /home/enomix/www/EsoterikoFinal/Maps/shps/nodes.shp  | psql -h localhost -p 5432 -d esoteriko -U postgres';
exec($command,$out,$ret);

Note: Don't forget to add the postgis extension to your postgresql database like:

CREATE EXTENSION postgis


texts: Transforming a image to a shape defined by image corners in earth coordinates - Python

texts: I have a source satellite image. It is a rectangular shape x by y. I also have the stored lon/lat coordinates of the corners and the center point in the xml meta file. Is there a way with the gdal or rasterio library in Python, to transform (scale, rotate, affine transform) the image to the shape defined in lon/lat coordinates of the corner points?


texts: Points proximity to polygons in ArcGIS for Desktop?

texts: I have a point feature class with 1725 records (spread out over NYC) and a polygon feature class with 13 records (also spread out over NYC). What tools or processes can I use to figure out what percentage of those points are within a 1/2 mile radius of any of the 13 polygons?


texts: Mosaic to New Raster outputs whole numbers instead of decimals

texts: I just ran the ArcMap Mosaic to New Raster tool to merge 5 rasters of pH levels in different farms. These individual rasters have values that are decimals. When I run the tool, the output is in whole numbers. Why is this occurring and how do I change it? Could it be an environment setting, a tool property setting, or a data property setting? Nothing looks even remotely out of place.


texts:
I create a map in QGIS, in which I added Grids with labels and Scale.
So in the map the scale is not equal to interline grids.
I know it is related to projection system and location. but how do I fix it in QGIS.

last I want to show labels of grids in Degrees not UTM. how?


texts: Error 99999 All raster processes fail in ArcGIS

texts:
I'm using ArcMap 10.8.1 and wanted to e.g. do a slope calculation on a Digital Elevation Model.
Now I´m getting Error 999999 which usually refers to a bug or some unexpected error. First I thought it´s the one DEM but no other DEM raster works. So the program has a problem here. Hillshade also doesn´t work, gives Error 010067, Aspect Error 99999, Contour lines work.
Raster calculator gives the error "format GRID" whatever that means now.
It seems like all map algebra processes are buggy for me.
Should I reinstall or do something else?

I´ve reinstalled ArcGIS on this computer: No Changes.
I´ve installed in on a second computer: No Changes.
The Elevation models used mainly for testing for the error are:
srtm_lev_bnd21.tif
Cell size: 0,00029550817, 0,00029550817
1.66 GB
32 bit
GCS_WGS_1984

SRTM_30.tif
Cell size: 30,30151988, 30,30151987
1.4 GB
16bit
WGS_1984_UTM_Zone_36N

jor_utm.img
Cell size: 29,23080851, 29,23080851
98.91MB
16bit
WGS_1984_UTM_Zone_36N

dem_l8_5x5_high_clip.img
Cell size: 0,5, 0,5
3.33 GB
32 bit
UTM_Zone_36N


texts: How to determine the centroid of polygons?

texts:
I have a polygon vector and I want to quickly attain the centroid of this polygon for insertion in a text document. I am working with QGIS. I searched through other questions and although some are somewhat close they don't quite answer my question.

Ideally, the centroid coordinates would be in the attribute table of the polygon and I could just copy and paste the coordinates.

I don't want to have to create a separate point vector file that represents the centroid and then find these coordinates and copy and paste them.


texts: Synchronizing data from QField to original SpatiaLite database

texts:
I'm trying to get started with QField. I created a SpatiaLite database with several SpatiaLite layers in my QGIS project. Exporting that project with QFieldSync is super easy, as is (offline) editing the data in QFIELD on my tablet. When I import and synchronize the data back to my QGIS project all the data is there and everything seems fine. However, the source path of the SpatiaLite layers has changed:
Original path: C:\Users\...\Documents\Demo.sqlite
Path after synchronization: C:\Users\...\QField\import\demo2\Demo.sqlite
The original data is thus not synchronized, instead a new file with the same name but in a different location is created.
My questions are:
Is this supposed to happen?
If yes, what is a smooth way to make sure the original data gets updated as well?
If no, what should I do differently?


texts: QGIS statistical summary gives exponential value

texts:
I need the precise sum of areas for a polygon layer, where I have the area values in an attribute. When I use the Statistical Summary window it shows the sum in exponential notation:
Statistic   Value
Count       232
Sum         6.68428e+06
Mean        28811.6
Median      19993.4

At first I thought it was a display issue related to column width, but no. And I hoped that copy-paste might reveal the true value, but also no.
The Area values have 11 decimal places. Can I get the sum to be the same degree of precision as the original data? (I don't actually need 11dp, but I would like to choose for myself the precision to report, rather than have QGIS dictate!)
There is an excellent answer to a previous related question  (Summing data from column in Attribute table) but it does not mention this issue. I am using QGIS 3.10.10.


texts: Clipping file geodatabase vector layer using QGIS

texts:
I have a Vector layer loaded from a .gdb (geodatabase directory). I'm seeing my contour lines fine.

I've drawn a polygon on top of the area I want to crop on a new
layer.
I select Geoprocessing Tools > Clip.
I input my .gdb layer as the Input Vector Layer and choose my new polygon layer, "block" as the clip layer. I save the file to the canvas and to the Output Shapefile, but nothing appears on the final drawing.

How might I go about cropping the feature class in this .gdb?

If I have to save it to a .shp file, how can I only convert the window view into the shapefile, rather than the entire layer?


texts:
I'm writing a function that returns routes using PgRouting (v.2.6)
The main problem I'm having is that I cant get the time less than 1 second per query.

The main query I need is get time, length, and the line, so I don't need all the returning rows but a sum() of them.

I have this query:

SELECT  sum(rutanodos.cost*60) as sum_cost_minutes, sum(routing.calles_2po_4pgr.km) as sum_km,
       round((st_length(st_union(st_transform(geom_way, 4258)))*100)::numeric, 3) as distancia_linea_km, st_asEWKT(st_union(geom_way)) as geom_json
   FROM pgr_dijkstra(
  'SELECT id as id, source source, target target, cost, reverse_cost FROM routing.es_2po_4pgr
   WHERE geom_way && ST_Expand(
        (SELECT ST_Collect(the_geom) FROM routing.es_2po_4pgr_vertices_pgr WHERE id IN (38785, 744540)),0.01)',
         38785, 744540, false, false) as rutanodos
             join routing.es_2po_4pgr on routing.es_2po_4pgr.id = id2

In which:

I obtain the route with pg_dijkstra() from a node to a node.
I join the nodes with the line.
Finally I group the list to get: SUM of costs, SUM of length, union of the geometry in EWKT.

And returns (in 1.1 second):

sum of cost minutes:345.276354
sum of Km: 643.367627099999
line length: 729.228 (to see the difference with returned km)
EWKT line geom: "SRID=4326;MULTILINESTRING((-1.7734135 42.2439494,-1.7732242 42.2408917),(-1.773393 42.2445657,-..."

After explain command:

- "Aggregate  (cost=8224.76..8226.67 rows=1 width=80)"
- "  ->  Nested Loop  (cost=0.68..8209.75 rows=1000 width=130)"
- "        ->  Function Scan on pgr_dijkstra rutanodos  (cost=0.25..10.25 rows=1000 width=12)"
- "        ->  Index Scan using calles_2po_4pgr_pk on calles_2po_4pgr  (cost=0.43..8.20 rows=1 width=126)"
- "              Index Cond: (id = rutanodos.id2)"

This works properly and with I get all the results that I need, but, as I said, it takes more than a second.

I've read some advices to optimize, but still can't get better time. For example.

I've created indexes in source, target, geom and a primary key.
Vacuum and cluster on index.
Created bounding box for passing only a part of the network instead of all the nodes.
My data comes from osm2po, with topology created.

I've read a lot of posts about, and I've seen that people can get times of 50ms, but the fastest time I can get is about a second.

Is there anything that I'm missing?

Anyone knows any other tip or recommendation to get better times?


texts: Seeking GIS based Customer Relationships Management (CRM) tools?

texts:
I'm searching for a tool like something around a "GIS based Customer Relationships Management Tool". The functionality I am looking for is:

collect: collecting customer related data in several formats like e-mails, postal letter as scans, documents, online-masks and so on..
store: store these files geographically referenced in a DB
categorize: attach several attributes over time to these customers like: wishes, needs, historical questions, business category ...
manage by query's: get reports according to queries via attributes or via geo-tag or via combinations like: "show all customers that like potatoes and mushrooms living in a buffer distance of 20km around point xy and related documents".

What are some useful tools meeting my needs to test?


texts: Converting Shapefile to GeoJSON - UTM6 and MTM3 versions give bad coordinates

texts:
I am trying to convert SHP to GeoJSON to work with Google Maps API's data layers.

I am working with ogre (I also tried mapshaper, but it always give bad coordinates).  When I convert SHP with WGS84 format, the conversion works fine.  However when I try to convert MTM3 and UTM6 SHP formats to GeoJSON, the coordinates are all screwed up.

Bad coordinates example: {"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[642014.75,4911335.859999999],[642005.96,4911335.92],[642005.6600000001,4911380.0600000005],[642003.6200000001,4911405.58] ...

I notice that on the ogre website, there are optional fields for Source SRS & Target SRS.  Am I suppose to add something there to get UTM/MTM conversions to work?


texts: Use csv file as raster file and stitch them together

texts:
I'm having an issue trying to stitch together 100 32bit Tiff files (coming from thermal imagery video and extracted as 32bit Tiff).

I used Hugin to mosaic the pictures extracted from the video which give good result but has the bad habits to distort the finale stripes.

I can as well extract as CSV, the file is then a matrix of 640 columns and 480 lines of temperature value. I was wondering if there is anyway to stitch those csv together in a way to obtain a no distorted finale file using R by exemple.

Is a script could be able to say ok I recognize this combination of value in those two csv and I can stitch them together.

We considered some structure from motion software to do this work but they usually don't handle the 32bit Tiff.

It is a bit unusual question but the development of UAV technology will lead to the multiplication of data as movie that will need to be transfer in GIS application so it might be interesting to find a solution.


texts:
I am trying to follow the solution posted at Converting projected coordinates to lat/lon using Python to convert NZTM coordinates into WGS84.
I wrote following script to covert NZTM to WGS84:
from pyproj import Proj, transform

inProj = Proj(init='epsg:')
outProj = Proj(init='epsg:4326')
x1,y1 = 1400677.021,4914233.543
x2,y2 = transform(inProj,outProj,x1,y1)
print x2,y2

I am not sure what would be the epsg number for NZTM. I tried to find it online but could not make it.
What is the EPSG number for NZTM?


texts: Precision and scale 0 after passing parameters

texts:
I have a field in GIS that it needs a precision and scale argument. Using Arcpy, I pass the argument like this:

d_fields = ["acres" , "p_dev"]
    for field in d_fields:
        arcpy.AddField_management('SE_for_DC.dbf', "temp","DOUBLE", 15,6 )
        arcpy.CalculateField_management('SE_for_DC.dbf',"temp","!{}!".format(field), "PYTHON_9.3")
        arcpy.DeleteField_management('SE_for_DC.dbf', field)
        arcpy.AddField_management('SE_for_DC.dbf', field,"DOUBLE", 15,6 )
        arcpy.CalculateField_management('SE_for_DC.dbf',field,"!temp!", "PYTHON_9.3")
        arcpy.DeleteField_management('SE_for_DC.dbf', "temp")

in where I create a temp field with the precision and scale parameters then copy the attribute over, then I recreate the new field with the right parameters and copy the data to it.

Whenever I open the dbf, both the precision and the scale fields have 0 values. How can I fix this error and why is it behaving like it?


texts: Checking field is empty and then adding new row with zeros?

texts:
I created a table with three fields but there are no values in these fields. Here is the screenshot of these fields.

I want to add a new row with zeros for all these three fields. Here is my code.

with arcpy.da.UpdateCursor(selectionname,("OBJECTID","FREQUENCY","SUM_Shape_Area")) as cursor:
for row in cursor:
   if row[0] == None:
      row[0] = 0
      cursor.updateRow(row)
   else:
      print selectionname + " is not empty"

Selectionname is my table, a parameter defined prior to this part of my code. First, I need to check if these three fields are empty, then adding a new row with zeros if they are empty. After several tries, I still could not add a new row to the table.


texts:
I´m trying to add a TAB-file containing text-objects to a SQLIte database. I´ve already added points, lines and polygons, but when I try to add texts I get an error - something about an erroneous action on the server (it´s in Danish, so I´m not quite sure how to translate it).

I´ve heard that I might be missing some additional columns or tables, but I cannot find any informations about it.

Any suggestions/comments are greatly appreciated.

Regards
Kaare


texts: Create Point based on location of minimum elevation inside a polygon in QGIS3

texts:
For QGIS3.0:
Using DEMs, I am drawing polygons around depressions. Inside that polygon, I need to generate not only the minimum elevation value but the location that value was extracted from; this location needs to be a point by the end.
The DEMs in use have near centimeter resolution over features that are a kilometer or so in diameter, therefore generating random points, or sampling all points methods are computationally intense; the PC is capable but future usage of this technique may be on lower-end computers.
Things tried:
Maximum altitude in polygons
How to plot a point on the highest elevation of a raster on a polygon. QGIS
Is there a way to find the minimum Z point feature inside of each polygon in a polygon layer?
The above-mentioned methods were modified from maximum they were written for into minimum, they return a value but not the precise location or require iteration over randomly derived samples of the feature.
I have been performing a similar technique in ArcGIS but trying to commit to QGIS once and for all. This is the last process in my conversion I have not been able to figure out.


texts: Importing symbology to specific unique value in ArcMap?

texts:
I am working with arcview 10.3. I trying to change the symbology of unique value "43" in layer "Total" to the same symbology as in layer "a2-land".

but in the layer properties when i choose import button i can't find option to do it

and as a result the whole "Total" layer's symbology turn to the same as layer "a2". The problem is that all other unique value in layer "Total" disappear and i do want that other values: 574,431,155, will remain the same as before:


texts: How to generalize polygon attributes in a regular polygon grid?

texts:
I'm trying to generalize the attributes of a set of polygons(let's say they're census blocks with a population property) into another, overlapping set of polygons(let's say it's just a grid with no properties)

I'm certain there's a word for this problem that I've never heard of and never thought to search... If that's so, please just point me in the right direction!

I want to make a choropleth map that shows the population density of the area without showing the shape of the census blocks themselves. I want to put their attributes into the regular grid. BUT, there's of course a lot of overlap, and irregular shapes don't fit neatly within a square.

I need something that can:

Assume that population is evenly distributed across the census block polygons.
for each square in the grid, calculate..

what polygons are even a little bit under it, and for each one, what portion of that irregular polygon is indeed under it.
Then, lets say we have an irregular census block polygon that has a population of 100, and is 37% covered by the grid-block of interest... It would need to calculate the amount of the population attribute to assign to the block, in this case, exactly 37, assign it to that block, and then do the same, adding for every other shape under that block before moving on.

Then do that for all the other polygons in the grid layer.

I should end up with the exact same total population, just abstracted and generalized into regular shapes.

How can I do this, particularly in QGIS?

I really just want to count up some point data (points in polygon) and create a choropleth corrected for population density.


texts: How to use a graphic for labeling in Qgis?

texts:
The "label" option in QGIS allows a wide range of possibilities for displaying attributes of an object in the map.

Is there a way to use a symbol/graphic (eg stopsign.png) for labeling objects of a (polygon) layer?


texts: Update from/to Z from attributes using arcpy Update Cursor

texts:
Trying to update Z values onto some point and line features from the level attributes on each feature.  My script works fine on points - the previous Z values of 0 now have the same values as the level field.  It doesn't work on my line features though.

I guess I can't write to the firstPoint.Z and lastPoint.Z on each line, however it is giving no errors at all (the script exits successfully).

pointFields = ['SHAPE@Z', 'Level']
lineFields = ['SHAPE@', 'FromLevel', 'ToLevel']

with arcpy.da.UpdateCursor(pointFC, pointFields) as cursor:
    for row in cursor:
        row[0] = row[1]
        cursor.updateRow(row)

with arcpy.da.UpdateCursor(lineFC, lineFields) as cursor:
    for row in cursor:
        row[0].firstPoint.Z = row[1]
        row[0].lastPoint.Z = row[2]
        cursor.updateRow(row)

Not all of my lines are 2-point lines, but I only want to update the first and last vertex of each line with the FromLevel and ToLevel values from the attribute table.  I'm trying to modify the code from How to update the Z value using an UpdateCursor? however I am trying to figure out how to update only the first and last points of my line.

Here is what I've attempted based on one of the answers from the other Q&A

with arcpy.da.UpdateCursor(lineFC, lineFields) as cursor:
    for row in cursor:
        newGeom = arcpy.Array()
        for part in row[0].getPart():
            newPart = arcpy.Array()
            for pnt in part:
                if pnt != None:
                    if pnt == row[0].firstPoint:
                        newPnt = arcpy.Point(pnt.X, pnt.Y, row[1])
                    elif pnt == row[0].lastPoint:
                        newPnt = arcpy.Point(pnt.X, pnt.Y, row[2])
                    else:
                        newPnt = arcpy.Point(pnt.X, pnt.Y, 0)
                    newPart.add(newPnt)
            newGeom.add(newPart)
        newShape = arcpy.Polyline(newGeom)
        row[0] = newShape
        cursor.updateRow(row)

But still it completes without error yet no Z values are present.

How can I update only the first and last vertices of my line with the Z values from attributes?

It does appear that if pnt == row[0].firstPoint: isn't working.  I did a print row[0].firstPoint and a print pnt and they return exactly the same values, yet it's not going into that if

print "First Point = {}".format(row[0].firstPoint)
print "Pnt = {}".format(pnt)

returns

First Point = 1875883.117 5826666.937 0 NaN
Pnt = 1875883.117 5826666.937 0 NaN

Yet it still hits the else.


texts: Getting Leaflet routing control result outside map in a div or table

texts:
All of the examples of Leaflet routing control display the result (iterenary) inside the map area.

I want to display it outside in a div or table element. There is ItineraryBuilder in the API but no example on how to use it.
My code
routing = L.Routing.control({
        waypoints: [
            L.latLng(37.56569, 22.8),
            L.latLng(37.56569, 22.7)
        ],
        routeWhileDragging: true,
        itineraryFormatter: "test"
    }).addTo(map);

I want the iterenary to be displayed in textarea below Path Data

    Path data
    <textarea class="form-control form-control-sm test" readonly id="path-results" rows="8"></textarea>


texts: Activate / Deactivate customized ITool from Form in ArcObjects?

texts:
I'm using ArcGIS 9.3/VB6, and I have a little issue with the activation of a custom Tool:

I have a Class, that implements both ICommand and ITool

When I click on ICommand, the form is launched
First, The CurrentTool is deactivated

Set Application.CurrentTool = Nothing

Then I need to activate the ITool, so I click on a button, I get my tool UID, and :

Dim pSelectedTool As ICommandItem
Dim pCommandBars As ICommandBars
Dim pApp As IApplication

Set pCommandBars = pApp.Document.CommandBars
Set pSelectedTool = pCommandBars.Find(pUID)
Set pApp.CurrentTool = pSelectedTool

The pApp.CurrentTool.Name contain the name of my ITool, but the problem is that the Tool is still deactivated, when I click on the Map nothing happen!

I found a question similar to mine, but with no help.


texts: QGIS Temporal support for WMS source showing latest image only

texts:
I'm trying to incorporate the sentinel definitive surface reflectance imagery into my project as sourced from the Digital Earth Australia OWS server (https://ows.dea.gov.au). I'd like to be able to view the different images over time using the Temporal manager, but QGIS seems to always be displaying the latest image rather than the image from the date range set in the Temporal manager.
QGIS version: 3.16
Data source: https://ows.dea.gov.au, product(layer): s2a_ard_granule_nbar_t
I have tried playing around with the temporal settings on the layer properties and then used the temporal manager to set the project temporal settings but no joy - it still just displays the latest imagery.
I have tried setting a static WMS-T temporal range to be just a single date (to force it to a particular date). That didn't work.
As an aside, I cannot seem to be able to set the "Pass project temporal range to provider" - it is telling me that the option is disabled because the project temporal range is not valid. But I've set the project temporal range to be a single date to keep things simple. I'm not sure if this is the heart of the problem.
I know the data source supports time based queries as I can access them via a browser. So I'm confident that the issue is in QGIS (either my configuration or an issue with the code).
The only way I've managed to hack this to work is to hard code a TIME=2017-08-25 parameter as part of the original definition of the WMS interface, and then tell it to "Ignore GetMap/GetFile URI reported in capabilities". This gets me the image from the specific date, but is not particularly workable if I want to be able to use the time slider in the Temporal manager to change the date.
Any advice on how to get this configured correctly?


texts:
Let me start by saying that I am not a GIS expert nor a cartographer ,  and my knowledge of GIS is limited to the quasi-null professional training I had as an architect.

My question here was triggered by another question I have asked on stackoverflow regarding the use of Google maps API.

My problem is this : I am making a web-app or a website, where I need to generate quasi random points on the surface of the earth.

Those points however , should be , like pointed on the SURFACE of the earth as intended in SOLID surface . (e.g. not in the middle of the ocean or deep inside a lake)
In shorts - I need to be able to "separate" the bodies of water from the solid surface .

I am a fairly okay programmer for web-technologies, especially PHP - so generating points (nodes) was a 1-minute job , but verifying that they are "accessible" points to humans without a floating vehicle was another issue.

I will be able to exclude coordinates from the generated array - but where can I get this info ?

Is such information (like worldwide bodies-of-water including seas oceans lakes) available as coordinates ? as boundaries ?

Does anyone here knows how to use the data from Google-maps for this purpose ? like I said in the original question - they obviously have this data (the bodies of water are blue) - but can users separate those layers ?

I have from previous works a fairly okay resolution (for my purpose)  of vector world map - But I have no Idea how to use it to generate even approximation coordinates ..

Please remember that I am not a GIS expert , so If you use any technical terms , Please try to explain them like you would explain to your 2-years-old neighbor´s kid .

Any help as to how to approach this will be more than appreciated.


texts: Is there a DOI or citation for the new release of SoilGrids OCS (released December 2019)?

texts: I would like to cite the OCS 0-30cm dataset from SoilGrids released in December 2019. Is there a DOI or pending publication that I can cite?


texts: How to set field values as labels in QGIS2 Print Composer

texts:
I am mapping the distribution of 325 bird species, displaying circles in the grid cells each species occurs in. The data contains as many duplicates of the NAME as there are cells in which it occurs. I produce a map by filtering for a NAME then want to print the maps using the NAME field in a label as a caption below the map. I would like to have it generated automatically rather than write it in manually - in anticipation of being able to generate all species maps as an Atlas (that's another issue).

So far, all I get by using [% name %] is [name], not the species name. Any ideas of what I might be missing? Nothing I've found by searching has helped so far.


texts: A simple way to merge lines separated by a space in QGIS

texts:
In QGIS after creating contour lines by extracting them from a DEM (Extraction > Contour) and smoothing them (Line Smoothing) I noticed this situation:

In order to smooth out the corners displayed in the screen more, I set the parameters of the tool in this way:

The result is this:

Every line is broken. Is there an easy way to join them without having to modify them one by one and without losing the data for all of them?


texts: Removing black box outside clipped raster layer using ArcGIS Desktop?

texts:
I have a problem.

I used the mosaic to raster tool to form one large raster of the study area.
Using the clip tool, I clipped the boundary of the study area to this large raster map.

BUT! There is now a black rectangular box outside the study area boundary that I cannot remove.

This black area is showing as 'nodata', so I cannot remove it using the symbology function.

How can I remove the black background using ArcGIS Desktop?


texts: Check validity of file geodatabase using arcpy

texts: AS the title suggests, how can I whether a GDB file is valid using arcpy? One approach is, only checking the whether the name ends with .gdb, but that's not a robust approach. Because a normal file can be ended with .gdb, where GDBs are essentially folders to OS. So how can I do that?


texts: ArcGIS - How to export only map borders to a transparent .png file?


texts: Printing map in Print Composer to specified name (not Print Composer 1, 2, etc)?

texts:
I am trying to print a pdf map, while using a .qpt template.

I can create a pdf output but I how do I name the created print composer?

Each time I run the tool it creates a new 'Composer 1' then 'Composer 2'...  and I need to reference the name of the composer in the print to pdf code for it to work (Composer 14 below).

This is my current code:

import os, sys
from qgis.core import*
from PyQt4.QtCore import*
from qgis.utils import*
from qgis.gui import*
import qgis
from PyQt4.QtXml import*

# Pan to selected feature
pan_feature = iface.mapCanvas().panToSelected()

#Next Set print composer to map canvas using python:
canvas = qgis.utils.iface.mapCanvas()
layers = canvas.layers()
myFile = '/Users/Michael/QGIS_template.qpt'
myTemplateFile = file(myFile, 'rt')
myTemplateContent = myTemplateFile.read()
myTemplateFile.close()
myDocument=QDomDocument()
myDocument.setContent(myTemplateContent, False)
comp = iface.createNewComposer()
comp.composition().loadFromTemplate(myDocument)

map_item = comp.composition().getComposerItemById('Map')
map_item.setMapCanvas(canvas)
map_item.zoomToExtent(canvas.extent())

# Next Print PDF using the printer plugin:
composertitle = 'Composer 14'
folder = '/Users/Michael/Capstone'
extension = '.pdf'
mp = qgis.utils.plugins ['MapsPrinter']
for composer in iface.activeComposers():
    title = composer.composerWindow().windowTitle()
    if title == composertitle:
        mp.exportCompo (composer, folder, title, extension)
        break

My plugin is useless if the composer Id need to be updated every time I run the tool.

These are the two main Q&As that I used to build my code above:

Set print composer to map canvas extent using python
Save Print/Map QGIS composer view as PNG/PDF using Python (without changing anything in visible layout)?


texts: Simplifying polygons without "trimming" them

texts:
I have a quite difficult issue: I need to simplify polygons, but without "trimming" them - no parts of polygon can be lost; furthermore, all holes in polygons should be filled (but this last is less important and can be done by hand). I presented an example of what I need - the brown polygon is original, and the red outline is desired extent.

There is a significant amount of similar questions on the StackExchange, but answers were not sufficient.

I use QGIS 2.18.3 and I made some trials using SimpliPy plugin - but results using tolerance: 10  are not so satisfactory. Second picture presents it: in some places work has been done well by the plugin, but some parts are excessively simplified. Also, there are parts that had not been simplified.

I wonder whether I should just try with SimpliPy settings and parameters - or should I look for another way (maybe GRASS / SAGA algorithms)


texts: Is it possible to download road network data from ArcGIS?

texts: I am trying to get the roads from ArcGIS for the Chicago region as a shapefile. I just want paths along the road. Is there anyway to export the information from ArcGIS, so I can use it?


texts: Python script from graphical modeller of QGIS 3.4 gives KeyError?

texts:
For several tasks in QGIS 3.4 I need a python script. As I am not used to use python I created the script with the graphical modeler (I modeled the tasks in the graphical modeler and exported the python script).
The first part of the script is basically

an object-filter to choose all points with a common attribute of a point layer
create a heatmap with all the points choosen with the object-filter

When I start the model from the graphical modeler it works as it should. But after exporting the python script and start this from the processing toolbox, it gaves a "KeyError", that the output of the object-filter is not specified.

The first two algorithms in the python script are as follow:

from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterVectorLayer
import processing

class Automatisierung_20190527(QgsProcessingAlgorithm):

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterVectorLayer('punkte', 'punkte', types=[QgsProcessing.TypeVectorPoint], defaultValue=None))

    def processAlgorithm(self, parameters, context, model_feedback):
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(8, model_feedback)
        results = {}
        outputs = {}

        # Objektfilter
        alg_params = {
            'INPUT': parameters['punkte'],
            'OUTPUT_gefiltert_1990': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['Objektfilter'] = processing.run('native:filter', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        feedback.setCurrentStep(1)
        if feedback.isCanceled():
            return {}

        # Heatmap (Kerndichtenschätzung)
        alg_params = {
            'DECAY': 0,
            'INPUT': outputs['Objektfilter']['OUTPUT_gefiltert_1990'],
            'KERNEL': 0,
            'OUTPUT_VALUE': 0,
            'PIXEL_SIZE': 100,
            'RADIUS': 15000,
            'RADIUS_FIELD': None,
            'WEIGHT_FIELD': None,
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['HeatmapKerndichtenschtzung'] = processing.run('qgis:heatmapkerneldensityestimation', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        feedback.setCurrentStep(2)
        if feedback.isCanceled():
            return {}

The error message says the following:

Results: {}
Traceback (most recent call last):
File "<string>", line 34, in processAlgorithm
KeyError: 'OUTPUT_gefiltert_1990'

Can somebody help me with fixing that error?


texts: Why are there blank in the ST_AsGeoJSON results?

texts:
I have a vector layer tipor several polylines.
I opened this layer in QGIS and exoportar to a spatial table in PostGIS.
ST_AsGeoJSON want to use the function to get the GeoJSON to draw on Google Maps.
I used the select below:

SELECT ST_AsGeoJSON (geom), Codlin
   FROM "lines";

I noticed that in the result column st_asgeojson few lines comes in white and records that return blank geom field of the table space is blank.
I think some records of my table are not georeferenced.
How can I fix this error?


texts: Removing lines that don't follow another

texts:
I have two line layers, the blue and the pink. I've already modified the pink one by snapping to the blue (i.e. there are pink lines that match the blue underneath), however now I want to get rid of the excess pink ones.
I've tried to remove them using v.clean, but using option 'rmdangle' just produces exactly the same layer.

I've also tried using a virtual layer with this query:

SELECT l1.*

FROM line AS l1, line AS l2

WHERE st_intersects(start_point(l1.geometry),end_point(l2.geometry))

OR st_intersects(start_point(l2.geometry),end_point(l1.geometry))

GROUP BY l1.id

but that gives me a weird result and definitely does not get rid of the excess pink lines.

Any suggestions?


texts: Label a feature based on Intersection with another feature class

texts: In Arcmap, I have a feature class that I want to label but only where there is NO intersection with another feature class in my MXD. I know one strategy would be to populate a field based on a spatial query with the other feature class first using geoprocessing but thats not really an option. Is it possible to do this in the SQL Label Query in Label Manager or a Python Expression? My data is stored in a GDB


texts:
I am new on spatial kernel density estimation with r and need some suggestions. Saying, I like to estimate the density for some event occurring at a location, for example, the probability of occurrence of a disease in each state, or the probability of soybean yield at each county. And then use the spatial kernel estimation result to calculate the probability that the occurrence will be lower than some specified level (threshold). What R package has the functions for spatial kernel density estimation?

The section 8.4.2 in the book, "Geographically Weighted Regression: The Analysis of Spatially Varying Relationships", have a topic about geographically weighted Kernels.This method seems to be what I need. Do anyone know any R package with functions for geographically weighted kernels?


texts: Coding python loop to select features of polygon layer than are completely contained in group of 11 layers using ArcPy?

texts:
This is my first attempt to code in Python.

I have one polygon layer, A, that covers the UK and another set/group of layers that each cover sections of the UK, B. I want to select by location features of A contained entirely within B. As B is split into lots of separate sections (due to its size) I need to do this in a loop.

So far I have this code which I adapted from a similar question of Automating Select by Location in ArcGIS Desktop?:

# Import system modules
import sys, string, os, arcgisscripting, arcpy

# Create the Geoprocessor object
gp = arcgisscripting.create()

mxd = arcpy.mapping.MapDocument("CURRENT")
df = mxd.activeDataFrame
groupName = "WBBS_waterways"

# Load required toolboxes...
gp.AddToolbox("c:/program files (x86)/arcgis/desktop10.3/ArcToolbox/Toolboxes/Data Management Tools.tbx")

# Script arguments...
A = "//FILES/UNIXHome/Git/waterbirds_project/data/OS_waterways.shp"
outfolder ="//FILES/UNIXHome/Git/waterbirds_project/data/"

basename = "OS_just_rivers"

allLayers = arcpy.mapping.ListLayers(mxd, groupName, df)[0]

for layer in(allLayers):
        print("A")
        D=layer.name
        # Process: Select Layer By Location...
        gp.SelectLayerByLocation_management(A, "COMPLETELY_WITHIN", D, "", "NEW_SELECTION")
        print("B")
        outname = outfolder + "\\" + basename + D
        print("C")
        #Process: Copy Features...
        gp.CopyFeatures_management(A, outname, "", "0", "0", "0")

And this is the error message I get:

Runtime error  Traceback (most recent call last):   File "<string>", line 28, in <module> ExecuteError: ERROR 000368: Invalid input data.

which I think means something is going wrong on this line: gp.SelectLayerByLocation_management(A, "COMPLETELY_WITHIN", D, "", "NEW_SELECTION"), it looks like python can't find layer D but I'm not sure how to make it recognise layer D.


texts: OpenLayers projection map issue CRS EPSG:32639

texts:
I would like to use CRS EPSG:32639 in map.

const initMap = (layerTypesList: LayerPackage[], id: number) => {
map = new ol.Map({
    target: 'map',
    layers: [...loadLayers(layerTypesList, id)],
    view: new ol.View({
        projection: 'EPSG:32639',
    }),
});
map.set('id', id);
const bounds: ol.Extent = transformExtent(GEODETIC_EXTENT);
map.getView().fit(bounds, { size: map.getSize() }); };

I have tried to use EPSG:4326, it works well however while using another CRS it gives errors as

  Unhandled Rejection (TypeError): Cannot read property 'G' of null


texts: ArcPy Graph function creates empty graph output?

texts:
Background:
I have some python code that I need help with parameter passing into the MakeGraph_management arcpy function. I would like to create a two-series vertical line graph for every pair of input point feature classes using the 'zip' function
(ie. input fgdb1/fc1 and fgdb2/fc1 >> MakeGraph >> graph1.png; input fgdb1/fc2 and fgdb2/fc2 >> MakeGraph >> graph2.png, ...).

My Problem:
Now I am trying to create graphs in batch using ArcPy/python 2.7.3, but have been struggling with the way these ArcPy functions work around the template file. I created a graph template in ArcGIS 10.2 with everything setup (axis labels, titles, footers, legend names, etc.) I just need to get the data input working to correspond with Series1 and Series2.

The code below successfully iterates and creates graphs for every pair of input datasets as planned. However, the graphs contain no plotted data. I thought I could just use MakeGraph_management and SaveGraph_management since I have a template instead of also working with the ArcGIS Graph class.

My Code:

# Import system modules
import arcpy
from arcpy import env

# enable overwrite outputs
arcpy.env.overwriteOutput = True

# input feature directory: yr0_transects
dir_yr0_transects = r"W:\Working\abcdefg_working\PointDistance\data\point_distance\year0\output\output_year0.gdb"

# input feature directory: yr1_transects
dir_yr1_transects = r"W:\Working\abcdefg_working\PointDistance\data\point_distance\year1\output\output_year1.gdb"

# output directory
outWorkspace = r"W:\Working\abcdefg_working\Profiles"

# Create list of year 0 Transect features
arcpy.env.workspace = dir_yr0_transects
yr0_transects_list = sorted(arcpy.ListFeatureClasses())

# Create list of year 1 Transect features
arcpy.env.workspace = dir_yr1_transects
yr1_transects_list = sorted(arcpy.ListFeatureClasses())

# Set workspace
arcpy.env.workspace = r"W:\Working\abcdefg_working\Profiles"

# Set local variables
graph_tee = r"W:\Working\abcdefg_working\Profiles\Profile_graph_template.tee"

for yr0_transect, yr1_transect in zip(yr0_transects_list, yr1_transects_list):
    try:
        #Cast to string because ListFC() returns unicode.
        string_yr0_transect = str(yr0_transect)
        string_yr1_transect = str(yr1_transect)

        # create title for graph
        graph_title = "Profiles, B.F.E." + string_yr0_transect[0:17] + string_yr0_transect[47:]

        # Create the graph object
        graph = arcpy.Graph()

        # Add a vertical bar series to the graph
        graph.addSeriesLineVertical(yr0_transect)
        graph.addSeriesLineVertical(yr1_transect)

        # Specify the title of the Graph
        graph.graphPropsGeneral.title = graph_title

        # Create graph output name
        output_graph_name = string_yr0_transect[0:17] + string_yr0_transect[47:] + "_profile" + ".png"

        # Output a vertical line graph, which is created in-memory
        arcpy.MakeGraph_management(graph_tee, graph, "temp_graph")

        #Execute SaveGraph
        mygraph = arcpy.SaveGraph_management("temp_graph", output_graph_name, "MAINTAIN_ASPECT_RATIO", "1024", "768")

        del mygraph

    except:
        print arcpy.GetMessages()


texts:
Can someone tell me what the projection is in this zone map from NOAA?

I am trying to replicate it (or get close) in QGIS.


texts:
How can I create a buffer with square ends along a line segment in qgis?

I've buffered lines using v.buffer in the GRASS plugin, but the -c option for square ends is not exposed in that interface.

Is there another way to do this with GRASS, or another tool?

Thanks!

update: I ended up using the grass qgis plugin with the "shell" module and just used the grass commands.  If you want to see the results in qgis you can use the "add grass vector layer" tool.


texts: Adding data from KML/shapefile/GeoJSON to OSM nodes

texts:
I have data of streets which is marked with the osm way/node ids. I would like to add speed limit data which is available in KML, shapefile, or GeoJSON to these OSM wayids (OSM wayid with start node and end node).

data.seattle.gov/Transportation/Seattle-Streets/jc8u-fewc

I was trying to convert the OSM map to KML (with ogr2ogr) but this failed in one way or another and now am experimenting with shapefiles.

How would you go about this problem?

I am new to map file formats I probably don't see the simple solutions.


texts: Renaming or deleting fields in multiple layers without creating new layers (QGIS)

texts:
I am very new to QGIS.

I have added three new fields using batch processor to all my layers (150 or so), but named them incorrectly. All I want to do is rename them correctly but I can't seem to do this. Everything involves creating new layers with the renamed column (refractor fields) or if I deleted the columns and add them again then again it creates new layers without them and keeps the old ones. It will take forever to save all these files under a new name and a new place, so I was wondering if there is any way to delete or rename the columns in the same layers and not create new ones.


texts: NDVI usage for ArcCollector and Crop Yield

texts:
So my question is as follows:

I am a GIS tech for a large vineyard management company. Every year a real struggle of ours is getting a crop estimation (with +- %15 accuracy, company wants it to be +-%5). I currently use 1)Historical Bunch Weights 2)Average amount of producing vines in a vineyard block 3) average Bunch numbers per vineyard block.

This has all been a very slow, painful experience with needing to send guys in the field with the Collector apps I create to manually count bunches, or producing plants vs dead/missing plants.

So my real question is that I am currently trying to use our NDVI's from last year (which I have in Arc already) to somehow extrapolate out the "producing vines" part of the equation. AKA I need to convert the NDVI to a line/point (?) feature that only contains producing vines (or dead/missing vines since that would be a smaller file). I then can send the field guys out to count bunches in pin-pointed regions of high vigor etc to estimate.

Any thoughts?


texts: Using Update Cursor to categorize crash ratios?

texts:
I'm new to coding and I need to somehow use an Update Cursor that I can convert the rates of crashes in a block group into 3 separate categories.  The crash rate is the difference between the number of crashes in a block group divided by sum of all road lengths in the block group.  The update cursor needs to go through my crash_rate field and then place them in my crash_category field.  Category 1 is low, category2 is medium, and category 3 is high.
My Max crash rate is 373.35998535, my average / mean crash rate is 51.8848051, and my lowest crash rate is 0.  Here's the code that I used to find the Max, Min, and Mean:

import arcpy
import sys
import os

fc = r"C:\GISAutomation\Project\GISAutoProject\GISAutoProject.gdb\F2F_BG_crash"

#field = "crash_ratio"
num_sum = 0
num_counter = 0
maximum = 0
minimum = 0
#cursor = arcpy.da.SearchCursor(fc,["crash_ratio"])
if arcpy.Exists(fc):
    print("exists")
with arcpy.da.SearchCursor(fc,["crash_ratio"]) as cursor:
    for row in cursor:
        num_sum = num_sum + row[0]
        num_counter = num_counter + 1
        if row[0]> maximum:  # doing the comparisons to find the maximum value
            maximum = row [0]
        if row[0]< minimum:
            minimum = row [0]
num_average = round(num_sum / num_counter,2)
print ("The average per blockgroup is", num_average, "miles")
print ("The minimum per blockgroup is", round(minimum,3), "miles" )
print ("The maximum per blockgroup is", round(maximum,2), "miles")
print ("script 1 finished")

Does anyone have any advice about how I can set my update cursor to place the different values into either high, medium, or low categories?  I'm really lost.


texts: Bounding box doesn't change the extent of the map

texts:
I want to display maps in Leaflet that are only visible on the bounding box. The code works fine without the bounding box but the maps are too many so displaying the maps is very slow. Here is my implementation of the solution that is not working.
var owsrootUrl = 'https://localhost.com/geoserver/somemap/ows';

var defaultParameters = {
    service: 'WFS',
    version: '2.0',
    request: 'GetFeature',
    typeName: 'Shambamap:Kenya',
    outputFormat: 'text/javascript',
    format_options: 'callback:getJson',
    SrsName: 'EPSG:4326',
    maxFeatures: 3000
};

var customParams = {
    bbox: map.getBounds().toBBoxString(),
};

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

var WFSLayer = null;
var ajax = $.ajax({
    url: URL,
    dataType: 'jsonp',
    jsonpCallback: 'getJson',
    success: function (response) {
        WFSLayer = L.geoJson(response, {
            style: function (feature) {
                return {
                    stroke: true,
                    fillColor: 'FF0000',
                    fillOpacity: 0.2,
                    opacity: 0.5
                };
            },
            onEachFeature: function (feature, layer) {
                if (feature.properties) {
                    var content = "
" +
                                    "
Plot No	" + feature.properties.plotno_1.toString() + "
" +
                                    "
FR Src	" + feature.properties.frsrcname + "
" +
                                    "
Locality	" + feature.properties.locality + "
" +
                                    "
Details
" +
                                    "
";
                    popupOptions = { maxWidth: 500 };
                    layer.bindPopup(content, popupOptions);
                }
            }
        }).addTo(map);
    }
});


texts:
I am reading some EVI values from 46 raster layers (250 m resolution) 231 rows, 481 columns each and applying Support Vector MEachine (SVM) to classify using the training data stored in excel file and finally writing back the predicted file into one raster.

The code is okay, but it's taking me around 100 hours to run in my Intel Core i7 (7500 U) 2.9 Ghz, 8 GB RAM machine

Can someone help me rewrite this code to optimize the performance?
The inner loop (j loop) in the following code is taking more than 80% of the total time.

library(xlsx)
library(e1071)
library(caTools)
library(gdata)
library(rgdal)
library(raster)
Sys.setenv(JAVA_HOME='C:/Program Files/Java/jre1.8.0_131')
library(rJava)
setwd("D:/PhD Dissertation/R/SVM/June15")

Read the tiff file with 46 bands

phenology<-stack("EVI_2005_All.tif")
names(phenology)
dimensions<-dim(phenology)
rowNumber<-as.integer(dimensions[1])
colNumber<- as.integer(dimensions[2])

Define new raster and its projection, extent and resolution

predict_raster<- raster(nrow=rowNumber, ncol=colNumber)
res(predict_raster)=250
predict_raster<- raster(nrow=231, ncol=481, xmn=163244, xmx=283494, ymn =3028497  , ymx = 3086247)
projection(predict_raster) <- "+proj=utm +zone=45 +datum=WGS84 +units=m"
print(predict_raster)
all <- read.xlsx("Training_data_2005.xlsx", sheetName = "EVI_2005_All")

Create a sample index of the total data set and divide it to train and test set

sample_index <- sample(466,350)
svm_train <- all[sample_index,]

When the class column is numeric in such cases:

Convert the class column to factor so that the SVM performs classification rather than regression

Then convert the training set into data frame

svm_train[,"Species"] <- factor(svm_train[,ncol(svm_train)])
svm_train <- as.data.frame(svm_train)
svm_test <- all[-sample_index,]

---------- Now Create an SVM Model with the training dataset--------------------

svm_model <- svm(Species ~ ., data = svm_train)
summary (svm_model)

p=0

for(i in 1:rowNumber)
{
  for(j in 1:colNumber)
  {
    #read the value for each 46 bands for particular coordinage i and j
    v <- (phenology[i,j])*0.0001

    #Now the data is ready to go for the model
    p <- as.integer(predict(svm_model,v))
    predict_raster[i,j]<-p

  }
  Sys.time()
  print("Working in Row #:")
  print(i+1)
}

print(predict_raster)
print(phenology)
writeRaster(predict_raster,"Predict_SVM_EVI_All_2005.tif","GTiff", overwrite=TRUE)
print("End time of  writing the output image:")
Sys.time()


texts:
I am new to GIS and am currently doing a project related to GIS.

Can somebody recommend me a good GIS software that can:

Create custom indoor maps
Do custom routes based on custom criteria
Generated map can be used for directions on mobile


texts: How to toggle layer visibility in ArcObjects?

texts:
ILayer.Visible seems to have no effect - How can I manually manage layer visibility and actually get the map control to update?

ILayer lyr = mapControl.get_Layer(i);

//this works for custom layers that I define but not for Bing_Maps.lyr or ESRI .lyr files that comes with ArcEngine or connections to ArcGIS Online.  For composite layers, I also tried iterating through sub layers and setting the visibility to false

lyr.Visible = false;
mapControl.Update();
mapControl.Refresh();
mapControl.ActiveView.Refresh();
mapControl.Invalidate();


texts: Selecting cells with values from DEM those are above specific value corresponding to each polygon?

texts:
i have a shapefile with 2000 polygons and a DEM. i want to calculate zonal mean for each polygon considering values from DEM which are above a particular value for each polygon.

i want to it with ArcGIS 10.2


texts:
I'have a raster image from Scran25 IGN. In Qgis, it opens perfectly but in R (with rgdal or raster packages), I do not know how to keep the original colors like this exemple?

I have tried

library('rgdal')
library('raster')

scan25<-readGDAL("./Banyuls/MyCreate/etudehydoAOC/etudehydoAOC/Scan25/vrt_4communes.vrt")
communes<-readOGR(dsn = "./Banyuls//MyCreate/Surface_4communes/",layer = "4communesAOC")

scan25<-raster(scan25)
test<-aggregate(scan25,10,mean)

plot(communes)
plot(test,add=T)

And the result is


texts: Definition query to display only part of layer that overlaps with specific polygon using ArcGIS Desktop

texts:
I have a very large roads data set and I would like to write a definition query that displays only the roads that are inside of a polygon. This would make my export to pdf a much smaller file size. I know I could clip them but I don't want to make a new dataset.

Also I could create layer from selection, but I am very curious if this possible with a definition query?


texts: How do I get the latitude and longitude of only the white pixels in an area?

texts:
I'm working on a project where we're trying to measure the shape of coastlines.  After some experimentation, I came up with an edge detector that has very little noise.  Now, my goal is to get the longitude and latitude of only the white pixels in an area and put them in a array or dictionary.  Then I can put them into MatLab or Excel and create a curve function.  Can some one help me with getting the long/lat or even an x/y value of only the white pixels?

Below is the code I have in GEE:

var test_curve =
/* color: #ffebcc */
/* shown: false */
/* displayProperties: [
  {
    "type": "rectangle"
  }
] */
ee.Geometry.Polygon(
    [[[17.936393079223645, -32.69828825968294],
      [17.936393079223645, -32.72363697964223],
      [17.970382031860364, -32.72363697964223],
      [17.970382031860364, -32.69828825968294]]], null, false);

// USGS Landsat 8 TOA Reflectance (Orthorectified) Set
var l8 = ee.ImageCollection('LANDSAT/LC8_L1T_TOA');

// Center on ROI and zoom to level 14
Map.centerObject(test_curve, 14);

/****************************************************
 * Preprocess the data set, filter by images related to
 * the ROI, within the 2017 year, and minimal cloud interferance
 * **************************************************/

// Specify only the area in the test_curve polygon
var spatialFiltered = l8.filterBounds(test_curve);

// Further filter the number of images by only
// choosing the ones from 2017
var temporalFiltered = spatialFiltered.filterDate('2017-01-01', '2017-12-31');

// This will sort from least to most cloudy.
var sorted = temporalFiltered.sort('CLOUD_COVER');

// Get the first (least cloudy) image.
var image = sorted.first();

/*******************************************************
 * Convert image to grayscale and normalize to reduce noise.
 * perform edge detection to isolate points on the curves.
 * *****************************************************/

// Create an NDWI image, define visualization parameters and display.
// This will reduce the noise and return and image with fewer artifacts.
var ndwi = image.normalizedDifference(['B3', 'B5']);
// this will reduce the palette to only black and white
// min: .50 returns the best result but they are hard to see
// so min: .41 has been used for visulization purposes.
var ndwiViz = {min: 0.41, max: 1, palette: ['000000','FFFFFF']};

// After experimentation, edge detection was done with laplacian8
// and magnitude set to .21, as it returned the cleanest results.
var p = ee.Kernel.laplacian8({ magnitude: 0.21, normalize: false});

// Apply edge-detection kernel
var edgy = ndwi.convolve(p);

/*******************************************************
 * Display image
 * ****************************************************/

Map.addLayer(edgy, ndwiViz, 'curve_1', true, 1);


texts: How do I set the GRASS_ADDON_BASE environment variable?

texts: I am using QGIS and now GRASS GIS but I am not overly computer savvy. I am trying to add a GRASS modules (r.lfp) to my GRASS setup. I am on a Windows computer and can use the GRASS shell (I used to use .DOS in the old days). I've tried installing the extension using g.extension: it installs successfully (to c:/users/Bill/.grass7/addons) but also says that it won't function until the GRASS_ADDON_BASE variable is set. Can someone explain how I do that in simple language? I also tried to install r.lfp by saying "r.lfp -s" but perhaps I don't have permission to do so since that didn't work either. I've also tried to add r.lfp through the wxpython gui but the computer just hangs up and doesn't fetch the extensions repository (does it need something called numpy?)


texts:
I have a standalone pyqgis program and it uses QgsMapCanvas as map display. I would like to set its extent to be empty. How do I do that?

I tried the following but it is not working, in the sense that the canvas's extent is not empty after setting an empty rect:

r = QgsRectangle()
if r.isEmpty():
   canvas.setExtent(r)


texts: How many function templates can be applied to single image service in ArcGIS for Server?

texts:
I am trying to minimize space used on our server by producing a single image service for each county in Maryland.  However, I need to display 6 different derivatives for each county: DEM in feet, meters, stretched, aspect, slope and hillshade.  Does anyone know how many function chain templates I can apply to a single image service? Or will I have to create a image service for each derivative?  The latter seems to defeat the purpose of being able to run on-the-fly processing to a mosaic dataset...

Right now I have each derivative as its own image service... Being that I am not even through half of the counties processing and I am out of memory, I need to find an alternative.


texts:
I have a lesson/class assignment that teaches the use of format and export environmental raster data, but the issue is the lesson is outdated from the new QGIS version.

Lesson states:

Multiply the values in two layers together, or (as here)
divide all values in a layer by 10. In QGIS, go to Raster > Raster Calculator. Now, select the layer you want to manipulate (bio_1) by double-clicking it, select the name and directory of the Output layer, and then type the expression to divide by 10 – use the buttons, not your keyboard to type / 10:

Then click OK!

The QGIS 3.0 program doesn't allow you to create a new output file (or to my knowledge).

What should I do?


texts: I get in QGIS: "ERROR[25]: Error reading raster. Possibly the raster path contains unicode characters"

texts:
I have done clipping, and classified the land according to the various characteristics in QGIS.
But when I try to save the processed image to my local machine I get the following error message:

"ERROR[25]: Error reading raster. Possibly the raster path contains unicode characters".

I don't know what to do.


texts: Extracting node id from OpenStreetMap data in QGIS?

texts:
There is a lines layer, and I want to extract node ids from the lines.

(the data can be downloaded from http://download.bbbike.org/osm/bbbike/Beijing/)

A possible answer is to extact vertices first (Vector -> Geometry Tools -> Extract vertices). However, it only generates some attributes, like vertex_index.

After that, I can only get its coordinates directly, but it is hard to retrieve the node id itself. (the osm_id in points).

To be specific, in the osm file, we can find:

<node id="5761117126" lat="39.9135365" lon="116.5852367" version="1"/>
<node id="5761117127" lat="39.9129296" lon="116.5880741" version="1"/>
<node id="5761117128" lat="39.9124111" lon="116.5877737" version="1"/>
<node id="5761117129" lat="39.9109294" lon="116.5873889" version="1"/>
<way id="4231222" version="1">
    <nd ref="25248785"/>
    ...
</way>

Now I need to get refs inside ways.


texts:
Although I've a regular r user I have only ever computed GWR in arcgis, and now I'm using R I'm just a little confused.

I seem to have got my GWR results and now need to go about plotting them on a map. I did this GWR on a spatial points data frame, which I would usually join spatially to a separate polygons. But with my results stored in an object called 'res.binomial', how can I then join them up to my other data?

EDIT: I've seen that results seem to be stored in res.binomial$SDF - so I guess the question is how I use those results??

bwG <- gwr.sel(pop~ demsat + lrscale + relig + euint + econ, data = newest, coords = locationsnew, verbose = FALSE)
bw.gwr.1 <- bw.gwr(pop~ demsat + lrscale + relig + euint + econ, data = newest, approach = "AICc",kernel = "bisquare", adaptive = TRUE)

DM <- gw.dist(dp.locat=locationsnew)

res.binomial<-ggwr.basic(pop~ demsat + lrscale + relig + euint + econ,data=newest, bw=bwG, dMat=DM, family ="binomial")

res.binomial


texts: Does Linux support "like '%c'" filter(cql_filter)

texts:
I can execute the following request:

http://My.host:8080/geoserver/wfs?request=GetFeature&service=WFS&version=1.1.0&typeName=topp:states&cql_filter=STATE_NAME like 'California'

But I can't execute the following one:

http://My.host:8080/geoserver/wfs?request=GetFeature&service=WFS&version=1.1.0&typeName=topp:states&cql_filter=STATE_NAME like 'C%'

Here goes the error message:
XML Parsing Error: no element found
Location: .http://My.host:8080/geoserver/wfs?request=GetFeature&version=1.1.0&typeName=topp:states&cql_filter=STATE_NAME%20like%20%27C%%27
Line Number 1, Column 1:

I think something is wrong when I use % he encode the request and something happens...

GeoServer:V2.1.4
Browser:FireFox
SO: Ubunto


texts: Why are some maps shown with polar coordinates (azimuthal) and some Cartesian (cylindrical)

texts:
I am very new here, and in no way am I a cartography expert like all of you here. I am currently doing a math paper on map projections, and I have noticed that some projections (mainly azimuthal) are mapped in polar coordinates when they could also be mapped in Cartesian, and vice versa goes for cylindrical projections.
I understand that this may make points easier to refer to, but are there any other more 'solid' reasons? I have researched many cartography related websites, and most of them just provide the mapping functions without justification, and some just mention that "... projection is most often shown in polar form".


texts: Identifying Lithuanian coordinate system?

texts:
I am trying to make sense of geodata which uses coordinates like the following:

(574457.93, 6058460.22)–(574382.72, 6055752.97)

These should be located just west of Vilnius, Lithuania. 54.642587, 25.151125 should be close to the middle between the two points. The points above should refer to the junctions enclosing that point.

The coordinates are labeled x and y, respectively. No further data (such as a Z coordinate or zone) is given.

I did not find any documentation which could give a clue about the coordinate system used. The Y component is close to the northing in meters I would expect, but I can’t make any sense of the X component. If it is easting in meters, the reference meridian would be around 9° E.

The data is from an open data portal of a Lithuanian government agency.

Any ideas?


texts: How to use the Intersects function in the Field Calculator of QGIS?

texts:
I would like to have an example on how to use the Intersects function in the Field Calculator of QGIS.

Basically, I want to know which feature from a polylines layer intersects a polygons layer and store the result in an attribute. (0 = not intersecting, 1 = Intersecting).

I'm aware of the geoprocessing tools "Clip" or "Intersect", but is it possible in the Field calculator ?


texts: Hausdorff distance (ST_HausdorffDistance) extremely slow for polygons?

texts:
I have a simple SQL statement for computing the Hausdorff distance between two polygons of the form:

SELECT ST_HausdorffDistance(geom1, geom2);

, which is taking a long time (for a pair of polygons with a little over 1000 vertices each). The full query (in WKT: https://pastebin.com/ndyxE0aD) takes almost 5 seconds on a computer with Intel i9 CPU (PostgreSQL 11/PostGIS 2.5 on GEOS 3.7.2).

I needed to do this computation many times to compute Hausdorff distances between all features from two polygon tables, each with about 100 features. (The pasted example are the first records in each table). The query took 3.5 hours for such a small dataset. Similar queries on much larger road tables (with 500 to 1000 features in each table) took less than a minute.

I was wondering why the ST_HausdorffDistance computation is slow on the polygon data?

Is it because of an implementation issue or bug of PostGIS or the GEOS function behind it? Or is there something wrong with the polygon data I use?

Has anyone had similar issues with the polygon Hausdorff distance?

-- Update --

I only used time command from shell to test the times. As pointed out in the comments, the actual CPU time isn't bad at all.

For one run, the real user and sys time are 4.6, 0.05, 0.01 seconds respectively. For a larger run involving two datasets of similar polygons (but with data directly read from database, not from ST_GeomFromEWKT), it's 138m8s, 0.2 s, 0.05 s, respectively. So the CPU time is indeed small.

The GEOS package is compiled from source:

CXXFLAGS=-std=c++17 ./configure --prefix /opt/geos --enable-python
make && sudo make uninstall && sudo rm -rf /opt/geos && sudo make install

My new question is:

Why is the real time (wall time) so long compared to CPU time? (I have a SSD.) and how to configure PostgreSQL to reduce the wall time?


texts: How do I get the building housenumber from OSM data?

texts:
I've made a map in which buildings in a city are displayed using shapefiles from Geofabrik and unfortunately, none of the shapes contain the housenumber.

When I download the .osm files from Geofabrik the multipolygon layer is empty, but when I download the OSM Data directly via the OSM download tool in QGIS, the multipolygon layer contains buildings, but only if the area where I download it is very small (roughly 30 buildings at a time).

Is there a way to get the multipolygons from OSM other then downloading 30 buildings at a time, then when those are downloaded, extract the housenumber from the other_tags field in the attribute table?


texts: Dissolve singlepart features into multipart based on another polygon?

texts:
I'm working on generating a water feature layer that eliminates a lot of the swamp/marsh areas in my county. Using raster analysis (iso unsupervised cluster, majority filter, boundary clean, etc) and converting to a polygon, I now have a layer that more accurately depicts surface water.

However, this also means that I now have a lot of spurious features. Here's a screenshot:

Green is my new layer and blue is the original. Those selected polygons, and all of the other green ones inside the blue boundary, should be 1 feature. However, the attribute table for this selection shows 25 different features. I want these small, isolated singlepart features to be merged together based on the original polygon that surrounds them. It would be like a spatial dissolve where you dissolve based on location instead of attributes. Unfortunately, the closest thing that I could find to my desired process was merging in the editor toolbar, and there's no spatial component to that--it's just selecting by clicking.

It seems like this should be achievable. Any suggestions?


texts: Filling/replacing nodata values of a raster layer in R

texts:
Is there a way to fill my NoData values in R?

I have a raster that has several nodata pixels, which looks like this:

(each white dot has no data)

For the moment, the no data values are assigned -Inf:

Formal class 'RasterLayer' [package "raster"] with 12 slots
  ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
  .. .. ..@ name        : chr "D:\\Users\cost_based_on_density_new.asc"
  .. .. ..@ datanotation: chr "FLT4S"
  .. .. ..@ byteorder   : chr "little"
  .. .. ..@ nodatavalue : num -Inf
  .. .. ..@ NAchanged   : logi FALSE
  .. .. ..@ nbands      : int 1
  .. .. ..@ bandorder   : chr "BIL"
  .. .. ..@ offset      : int 0
  .. .. ..@ toptobottom : logi TRUE
  .. .. ..@ blockrows   : int 1
  .. .. ..@ blockcols   : int 892
  .. .. ..@ driver      : chr "gdal"
  .. .. ..@ open        : logi FALSE
  ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
  .. .. ..@ values    : logi(0)
  .. .. ..@ offset    : num 0
  .. .. ..@ gain      : num 1
  .. .. ..@ inmemory  : logi FALSE
  .. .. ..@ fromdisk  : logi TRUE
  .. .. ..@ isfactor  : logi FALSE
  .. .. ..@ attributes: list()
  .. .. ..@ haveminmax: logi FALSE
  .. .. ..@ min       : num Inf
  .. .. ..@ max       : num -Inf
  .. .. ..@ band      : int 1
  .. .. ..@ unit      : chr ""
  .. .. ..@ names     : chr "cost_based_on_density_new"
  ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
  .. .. ..@ type      : chr(0)
  .. .. ..@ values    : logi(0)
  .. .. ..@ color     : logi(0)
  .. .. ..@ names     : logi(0)
  .. .. ..@ colortable: logi(0)
  ..@ title   : chr(0)
  ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
  .. .. ..@ xmin: num -15.1
  .. .. ..@ xmax: num -7.64
  .. .. ..@ ymin: num 7.19
  .. .. ..@ ymax: num 12.7
  ..@ rotated : logi FALSE
  ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
  .. .. ..@ geotrans: num(0)
  .. .. ..@ transfun:function ()
  ..@ ncols   : int 892
  ..@ nrows   : int 658
  ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr NA
  ..@ history : list()
  ..@ z       : list()

How can I assign a new value (ex: 10000) to these nodata pixels in R?


texts:
I am generating geographical services from a PostgreSQL/PostGIS database in GeoServer using python, more specifically the gsconfig library.

I am using commands such these demonstrated in this link: https://www.gaiaresources.com.au/geoserver-rest-python/

So far so good. The problem is that when I generate a geographical service I am creating both WMS and WFS. And I just want, in some cases, to generate a WMS but not a WFS.

How can I do this configuration?


texts: Reproject Layer without losing Layer Style

texts:
I'm trying to reproject several layers via batch process from one CRS to another. In QGIS 3.6 it seems to work but when I reproject all my layers the layer style gets lost. Every layer has a different style.

Is there any way to reproject the layers without losing the layer style?


texts: How to add style in GeoServer with restsharp

texts:
I can add new style application/zip with Postman into GeoServer,  but I can't do it with  C#-RestSharp.
How can I do this?
      string url = "http://" + server + "/geoserver/rest/styles?name="+ StyleName;

            var client = new RestClient(url);
            var request = new RestRequest(Method.POST);
            string authInfo = username + ":" + password;
            request.AddHeader("authorization", "Basic " + Convert.ToBase64String(Encoding.ASCII.GetBytes(authInfo)));
            request.AddHeader("cache-control", "no-cache");
            request.AddHeader("content-type", "application/zip");
            request.AddHeader("authorization", "Basic YWRtaW46Z2Vvc2VydmVy");
            FileStream uploadFileStream = System.IO.File.OpenRead(@"D:\aa.zip");
            request.AddParameter("application /zip", uploadFileStream, ParameterType.RequestBody);

            IRestResponse response = client.Execute(request);

   string url = "http://..../geoserver/rest/styles?name=" + StyleName;
   wsn = "mystyle";

   WebRequest request = WebRequest.Create(url);

   request.ContentType = "application/vnd.ogc.sld+xml";
   request.Method = "POST";
   string authInfo = username + ":" + password;
   request.Headers["Authorization"] = "Basic " +
   Convert.ToBase64String(Encoding.ASCII.GetBytes(authInfo));
   request.Credentials = new NetworkCredential(username, password);
   byte[] buffer = Encoding.GetEncoding("UTF-8").GetBytes(wsn);
   Stream reqstr = request.GetRequestStream();
   reqstr.Write(buffer, 0, buffer.Length);
   reqstr.Close();


texts: Which type of point cloud is more accurate: from airborne LiDAR or aerial photos?

texts:
I have been offered two types of point cloud for my 3D Project. The data will be used by the environmental agency in various projects e.g. flood-risk analysis, dyke conservation, and nature conservation.

The first product is from airborne laser scanning and the second from aerial photogrammetry. The accuracy of the data is more important, meaning that +-20cm rounded to mm is more important than having +-1m with 6 decimals places.

Which of these two products would be more accurate in the sense of "correct"?


texts: Updating PostgreSQL geometry column using SQL (Trgger)

texts:
I'm looking to inderstand the role of date_part('epoch'::text, now()) WHERE id=NEW.id; in the trigger function below.

CREATE OR REPLACE FUNCTION fn_test_table_geo_update_event() RETURNS trigger
AS $fn_test_table_geo_update_event$
  BEGIN
    UPDATE test_table SET
    geom = ST_SetSRID(ST_MakePoint(NEW.longitude,NEW.latitude), 4326),
    updated_ts = date_part('epoch'::text, now()) WHERE id=NEW.id;

    RAISE NOTICE 'UPDATING geo data for %, [%,%]' , NEW.id, NEW.latitude, NEW.longitude;
    RETURN NULL; -- result is ignored since this is an AFTER trigger
  END;
$fn_test_table_geo_update_event$ LANGUAGE plpgsql;


texts: How to copy layers between MXD's with arcpy?

texts:
In ArcGIS I need to being able to have a couple of MXD documents open and copy-paste groups of layers from one MXD's table of contents to another with python script.

Is this possible and, if so, how?


texts: Training manual 6.2.9: Feature (1) has invalid geometry. How do I troubleshoot this?

texts:
I've stumbled upon the following error referring to invalid geometry in Feature 1. What is feature 1 and how do I troubleshoot this error?
Loading resulting layers
Algorithm 'Extract by location' finished
QGIS version: 3.12.3-București
QGIS code revision: 8234261527
Qt version: 5.11.2
GDAL version: 3.0.4
GEOS version: 3.8.1-CAPI-1.13.3
PROJ version: Rel. 6.3.2, May 1st, 2020
Processing algorithm…
Algorithm 'Extract by location' starting…
Input parameters:
{ 'INPUT' : 'C:/Users/JOEBLOGS/Documents/vector_analysis.gpkg|layername=Buildings_34S', 'INTERSECT' : 'C:\\Users\\JOEBLOGS\\Documents\\vector_analysis.gpkg|layername=road_schools_buffer_intersection', 'OUTPUT' : 'ogr:dbname=\'C:/Users/JOEBLOGS/Documents/vector_analysis.gpkg\' table=\"well_located_houses\" (geom)', 'PREDICATE' : [0] }

Feature (1) has invalid geometry. Please fix the geometry or change the Processing setting to the "Ignore invalid input features" option.
Execution failed after 0.28 seconds

Loading resulting layers
Algorithm 'Extract by location' finished


texts: Extend Lines To Another Feature (Polyline)

texts:
I am trying to convert the polyline shapefile shown in the picture below to polygon in order to fill in the "narrower" area between them as I am trying to get the right-of-way shapefile.

However, two problems arise:

The Feature to Polygon will automatically fill the bigger space between the lines.
The lines aren't exactly connected. The gaps are shown below

These lines are divided by different tiles. Each polyline shapefile is in a grid shown as below.

What I want to do is to first extend the lines to the grid boundary as there is also gap between them. If that works, then run a iterator to all the polyline features. Then merge them all.

Is there a way for me to extend the lines to the grid boundary?


texts:
I'm working with a feature class where I've generated 11 match tables for attachments for 11 chronological and numerically sorted fields.  After going through the attachment process, everything is sorted correctly in the feature class, but when it's published to our ArcGIS Online webmap, it's completely unsorted in the attachment section of the pop-up box.

Is there a reason for this or controls for the pop-up box to sort the attachments numerically?

I've attached screenshots of the same polygon's attachments. The sorted list is off the SDE, the unsorted list is the published pop-up box from our webmap.


texts: VectorLayer still empty even if attribute has been set with value

texts:
I am working with the QGIS API C++ 2.8.2, my application is written in QT 4.5.3 / C++

I had to create a VectorLayer in memory with attributes and features and export it to a shapefile to import it in QGIS application

So the creation of the VectorLayer works, it is based on the QGIS Python programming cookbook page 56, then  I decide to add attributes and set features.

Here is the code for the attributes:

QgsVectorDataProvider * lDataProvider = lVectorLayer->dataProvider();
QDateTime MyTime;
lVectorLayer->startEditing();
// add attribute
QList <QgsField> lMyAttribField;
lMyAttribField << QgsField("Id",QVariant::Int)
<< QgsField("Label",QVariant::String)
<< QgsField("Type",QVariant::String)
<< QgsField("Color",QVariant::String)
<< QgsField("PosX",QVariant::Double)
<< QgsField("PosY",QVariant::Double)
<< QgsField("TimeStamp",QVariant::String);

lDataProvider->addAttributes(lMyAttribField);

It is based also on stackexchange sample and cookbook QGIS
Then I set the feature like this:

// add feature
QgsFeature MyFeature;
QgsFeatureList MyFeatureList;
MyFeature.initAttributes(7);
MyFeature.setGeometry( QgsGeometry::fromPoint(QgsPoint(QPoint(5.9,43.12))));

for (int i=0; i<2; i++) {
    QgsAttributes MyAttributes;
    MyAttributes << QVariant(0) << QVariant("Test") << QVariant("Neutre") << QVariant("#556677")
    << QVariant(5) << QVariant(43) << QVariant(QString::number(MyTime.toMSecsSinceEpoch()));
    // set the feature with the attributes value
    MyFeature.setAttributes(MyAttributes);
    MyFeature.setFeatureId(i);
    MyFeature.setValid(true);
    // add feature in the list
    MyFeatureList << MyFeature;
}
// add the feature to the dataprovider
lDataProvider->addFeatures(MyFeatureList);

then finalise the edition of the VectorLayer and save to a file:

// sauvegarde des  attributs
lVectorLayer->commitChanges();
lVectorLayer->updateExtents();
lVectorLayer->updateFields();
// add to the registry
QgsMapLayerRegistry::instance()->addMapLayer(lVectorLayer);

QgsVectorFileWriter * writer;
writer = new QgsVectorFileWriter("/tmp/myshapes.shp",
    "CP1250",lVectorLayer->pendingFields(),
    QGis::WKBPoint,mpMapCanvas->mQgsCoordinate,
    "ESRI Shapefile");

    qDebug() << "Count of feature =" << lVectorLayer->featureCount();
    // insert the vector layer into the vectorlayermap
    mMapLayer.insert(pType, lVectorLayer);

When I launch the  application, I load a S57 file, map has been correctly displayed in my canvas ,I have a button to add Track or Layer, when I click on it everything seems to works, file .shp has been created in /tmp.

So I launch QGIS application, I choose Add vector, select my shapefile, the file appears without error in the layer list, when I open the attribute table there is only one field!!!! and no features!!!!!


texts: Im using QGIS 2.6.1 Brighton on an IMac with OSX 10.8.5. QGIS 1.8 worked on the same machine, as did QGIS 2.0. AFter installing 2.6.0 and then 2.6.1 I opened a file I'd used before and several queries were done. The attribute table had the selected items, but they did not appear on the screen as they had in the earlier version. When I tried to add the saved file to the map, it did it once and never again. The error message was that the shp file was not a valid data source. The shapefile does not appear to be saved. What am I missing?


texts: how to finds draw extents includes my point extent in openlayers3

texts:
i draw a DragBox on openlayers map and when i use drawend event, i can see thats extent.
and i have some point on my map. my point is selected in EXTENT want to check whether with using openlayers3 api

draw = new ol.interaction.DragBox({
        /* dragbox interaction is active only if alt key is pressed */
        condition: ol.events.condition.altKeyOnly,
        /* style the box */
        style: new ol.style.Style({
            stroke: new ol.style.Stroke({
                color: [0, 0, 255, 1]
            })
        })
    });
    /* add the DragBox interaction to the map */
    map.addInteraction(draw);

    draw.on('boxend', function (evt) {

       //if(evt.target.getGeometry().getExtent()**.isContains(MYPOINT.extent)**)

    });


texts: How to get around "The field is not nullable" error?

texts:
I have a shapefile of lines.
I want the angle of every line calculated in a field.

I am using the script below in arcmap, but when it encounters a line that is horizontal, it stops calculating any remaining values and gives the following error:

  The field is not nullable.

I searched for answers to similar questions, not getting too far with them.

Where I originally found this code:
https://geonet.esri.com/thread/20473#21719

Pre-logic Script code:

import math
def GetAzimuthPolyline(shape):
 radian = math.atan((shape.lastpoint.x - shape.firstpoint.x)/(shape.lastpoint.y - shape.firstpoint.y))
 degrees = radian * 180 / math.pi
 return degrees

Expression:

GetAzimuthPolyline( !Shape!)


texts: Why do I need to restart tomcat when using WFS from GeoServer across different apps and workspaces?

texts:
I've built a couple of apps that consume WFS 1.1.0 from GeoServer 2.2.2. I organise my layers in different workspaces, and they have a range of sources including shapefile, oracle, postgres etc.

I've noticed that wfs services only work in one application at a time. If i bounce tomcat then whichever app I go to first will be able to consume wfs. If i hit a second app, no luck. I can see the request for the second app being made in the logs:

14 May 16:06:24 INFO [geoserver.wfs] -
Request: getFeature
service = WFS
version = 1.1.0
baseUrl = http://localhost:8080/geoserver/
query[0]:
    srsName = EPSG:3857
    typeName[0] = {http://www.openplans.org/usaid}jamoats_geom
outputFormat = text/xml; subtype=gml/3.1.1
resultType = results

Puzzled. Am i missing some config or security setting?

Thanks!

<---Update after iant's request for more info-->

I inherited this stack from another developer...
I'm running nginx/0.7.65 & Tomcat 6 on Ubuntu 10.04.4 LTS.
geoserver.war is deployed at /var/lib/tomcat6/webapps.
I build a simple web-app (js/css/html) and package as a .war.
I psftp this over to my server next to geoserver.
I add sites to nginx for each web-app by editing /etc/nginx/sites-available/default.

I interrogated the XHR responses for a successful and an unsuccessful WFS request. I'm using two projects here to debug; 1.climvar and 2.usaid. If I load up climvar after restarting Tomcat my wfs features display. If I load up climvar after restarting Tomcat, BUT visit usaid first, no wfs features display. I looked at the responses in more detail. In both responses 82 features came back:

In the successful response:

xmlns:climvar=”http://www.fasga.com/climvar"...

In the unsuccessful response:

xmlns:usaid=”http://www.fasga.com/usaid”.

So it seems the wrong namespace is being 'remembered' across my WFS. There is no difference between the request payloads. In the unsuccessful response, all the prefixes are null. Instead of:

<climvar:id>0</climvar:id>

I get:

<null:id>0</null:id>

Here's my js layer definition for climvar:

    var layerOpg = new OpenLayers.Layer.Vector(nsClimVar.Constants.OPG_LAYER_NAME, {
        //id : "layerOpg",
        styleMap : nsClimVar.Styles.StyleMaps.styleMapOpg,
        strategies : [new OpenLayers.Strategy.Fixed()],
        protocol : new OpenLayers.Protocol.WFS({
            version : "1.1.0",
            url : "http://" + nsClimVar.Constants.SERVER_IP + "/geoserver/climvar/wfs", // geoserver wfs endpoint
            featurePrefix : "climvar", // geoserver worspace name
            featureType : "opg_new", // geoserver Layer Name
            featureNS : "http://www.fasga.com/climvar",
            srsName : "EPSG:900913" // default
        })
    });

Layers in usaid are declared in the same way. In GeoServer, both  workspaces have the namespaces set correctly, both are enabled, both have WFS enabled. I've tried removing the featureNS and featurePrefix params but this makes no difference.

The only workaround I can think of is having all my layers in one workspace which I don't want to do.

Any ideas?


texts: How to select features where attribute is not NULL?

texts:
I have an OS AddressBase shape file. One of the attributes is 'paoText'. It is either a text string like 'Anchor House' or is NULL.

I want to extract all the records where there is a text string i.e. the attribute is not NULL.

I have tried Advance search with paoText = NULL with the idea of inverting the result to get all those records that contain text i.e. a property name, but my search returns nill records.

I know there are records with NULL in that attribute because I can see them. What am I doing wrong please?


texts: Debugging AttributeError: ResultObject: Get attribute save not found

texts:
Below is the script for batch clipping of multiple raster images using a single shapefile.

import arcpy, glob, os
from arcpy import env
from arcpy.sa import *

# Set the input workspace
arcpy.env.workspace = r"F:\madhavi\images_to_clip"

arcpy.CheckOutExtension("Spatial")

# Absolute path to your mask layer
mask = r"F:\madhavi\shapefile\shp_gang_only_final.shp"

# Copying all the input rasters in an array
rasters = arcpy.ListRasters()

# Loop through rasters, append names and save files
for raster in rasters:
    output_raster = raster.replace(".TIF", "_clip.TIF")
    rasterObject = arcpy.gp.ExtractByMask_sa(raster, mask, output_raster)
    rasterObject.save(r"F:\madhavi\clipped_images")

There are two problems that I am encountering with this script.

Firstly, the first output image is getting saved in the folder where the input images are stored and not in the desired folder where I want to save the output images.

Secondly, the for loop is not working beyond the first input image.

The screenshot of the error is as follows:


texts: I am trying to use growler into a plugin, does anybody know how to set up a working growler or can point to some documentation?