texts: Creating DEM from point data in CSV files using ArcGIS Desktop?

texts:
I have a lot of csv files which contain point data containing height, aspect and slope measurements at a resolution of 12,5 meters.

I need to create a dem from this data and have ArcGIS 10 with spatial analyst and 3d Analyst. I did this before with 9.3 but I am not sure if my methodology is correct.

Create point shapefiles from the csv file
interpolate to raster (here I never know whether to use spline, kriging or some other method)?

Last time I did something wrong and rounded all the values up to 1m even though the data was to cm accuracy so this time I want to do it right!


texts: Concatenate string field in spatial join in QGIS

texts:
I'm doing a spatial join between a point layer and a polygon layer. The point layer has a field called "type" which has string values. When I join this to the polygon layer, I would like to get a field concatenating all the "types", in addition to the count of points in each polygon.

Here's an image of the point attribute table:

And here's an image of the two layers:

I know this is possible in ArcGIS with "Output field properties" and "Merge rules", but I can't figure out how to do this in QGIS.


texts: How to polygonize raster with range of values?

texts:
I have a raster with range 0 - 1. I want to polygonize them in QGIS according to their value to be later classified from "Very Low" to "Very High".

raster_value
possible_polygon_value
class

0-0.2
0
Very Low

0.2-0.4
1
Low

0.4-0.6
2
Medium

0.6-0.8
3
High

0.8-0.1
4
Very High

I am guessing there is a way to do it by using Raster calculator first to assign values to each pixel first. However, how? There is no "IF" and "ELSE" statement in raster calculator.
What I think:
IF ("raster@1" <= 0.2) * 0,
ELSE IF ("raster@1" > 0.2 AND "raster@1" <= 0.4) * 1,
ELSE IF ("raster@1" > 0.4 AND "raster@1" <= 0.6) * 2,
ELSE IF ("raster@1" > 0.6 AND "raster@1" <= 0.8) * 3,
ELSE IF ("raster@1" > 0.8) * 4


texts: How to create new service from existing service

texts:
I have an existing ArcGIS 10.0 Server, and have just setup an ArcGIS 10.4 Server. How do I go about moving a particular service from 1.0 to 10.4. I can locate the related mxd file as well as the SQL SDE database, but aside from that I don't know where any other data files are.

I am a developer who has been elected to do this admin side of GIS without any prior knowledge.

I really don't know where to begin aside from calling up ESRI and saying "help me".


texts: QGIS/ArcGIS projection problem between ESPG:4616 and ESPG:30165

texts:
I have two layers, one in ESPG:4616 (JGD2000) and one in ESPG:30165 (Tokyo, Japan Plane Rectangular CS V). In QGIS, layer are not matching as you can see below.

I tried everything in QGIS, modifying the project projection, with and without "on the fly", recording the layers in another coordinate system, etc. At this point, I through that the problem came from the data.

However, when I open both the layers in ArcGIS, the software do an automatic transformation (one projection into one another) and layers are matching as you can see below.

At this point I am lost, I even tried to record both the layers in a brand new projection in ArcGIS, still layers are not matching in QGIS.

Here is a link with a small sample of both layers

Dropbox file

UPDATE:

ArcGIS also has a problem, when I add a basemap, we can see that both layers are not overlapping

While using open layer in QGIS, one layer is overlapping the background map.


texts: QGIS 3 - Python : QgsVectorLayerJoinInfo issue

texts:
I'm trying to use Python to join a table to a vector layer within QGIS (v. 3.2.0)

Here is the script I use :

>>>root = QgsProject.instance().layerTreeRoot()
>>>ids = root.findLayerIds()
>>>target=QgsProject.instance().mapLayer(ids[1])
>>>print ("target :",target.id())
target : du_fl_pts_e154c28b_4812_41a2_9a78_bf95b7a97926
>>>layerToJoin=QgsProject.instance().mapLayer(ids[0])
>>>print("Layer to join :",layerToJoin.id())
Layer to join : TAXREF_9ab461d0_81d1_4095_b1d0_dce90b20600e
>>>lien = QgsVectorLayerJoinInfo()
>>>lien.targetFieldName = 'CD_REF'
>>>lien.joinLayerId = ids[0]
>>>lien.joinFieldName = 'CD_NOM'
>>>target.addJoin(lien)
True

But after the execution of this script, I have no join on my target layer.

Interesting fact : when I test the same script on the same layers/tables on QGIS 2.18.21 (just changing QgsProject -> QgsMapLayerRegistry and QgsVectorLayerJoinInfo -> QgsVectorJoinInfo) the join works.

Any idea why this doesn't work on QGIS 3?


texts:
I'm trying to create an algorithm similar to minimum bounding box (though it may end up looking nothing like it). In this case the angle will be passed as a parameter and given the angle I need the smallest rectangle covering all my points/polygons. So far my line of thought is to find the center of my points (centroid algorithm), and from there create two parallel lines with same angle as paramter angle, and two more lines perpendicular to them. Then using iteration move these lines outwards (in opposite directions) until they contain all points. Also doesn't have to be exact minimum bounding box, approximate works (I guess would depend on size of each iteration step).

Here is my code so far. I have already dissolved all my polygons into one. I then take a convex hull to reduce vertices. I then put all vertices in a list - not sure if this helps yet...

a = layer.getFeatures()
for feat in a:
    geom = feat.geometry()
a = geom.convexHull()
vertexId = QgsVertexId()
vertices = []
b = a.constGet().nextVertex(vertexId)
while b[0]:
    vertices.append(b[1])
    b = a.constGet().nextVertex(vertexId)

Notes: At some point I need to pass the angle of box. I am using QGIS 3, and need to create this in Python. Layer 'layer' has one geometry, the dissolved polygon of all other polygons - maybe iteration isn't needed to access it..

Please let me know if I should be passing on more details/info.


texts: geopandas read_file breaks in fiona open()

texts:
I am trying to read-in a shp file and can't seem to do it with the geopandas version 0.2.1. Here is the error:

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-28-b4165d16b7f3> in <module>()
----> 1 gpd.GeoDataFrame.from_file('D:/NHDPlusV21/NHDPlusMS/NHDPlus06/NHDSnapshot/Hydrography/NHDFlowline.shp')

C:\Users\Rdebbout\AppData\Local\Continuum\Anaconda2\envs\nusc\lib\site-packages\geopandas\geodataframe.pyc in from_file(cls, filename, **kwargs)
    173
    174         """
--> 175         return geopandas.io.file.read_file(filename, **kwargs)
    176
    177     @classmethod

C:\Users\Rdebbout\AppData\Local\Continuum\Anaconda2\envs\nusc\lib\site-packages\geopandas\io\file.pyc in read_file(filename, **kwargs)
     26         else:
     27             f_filt = f
---> 28         gdf = GeoDataFrame.from_features(f, crs=crs)
     29
     30     return gdf

C:\Users\Rdebbout\AppData\Local\Continuum\Anaconda2\envs\nusc\lib\site-packages\geopandas\geodataframe.pyc in from_features(cls, features, crs)
    185         """
    186         rows = []
--> 187         for f in features:
    188             if hasattr(f, "__geo_interface__"):
    189                 f = f.__geo_interface__

fiona\ogrext.pyx in fiona.ogrext.Iterator.__next__ (fiona/ogrext.c:20633)()

fiona\ogrext.pyx in fiona.ogrext.FeatureBuilder.build (fiona/ogrext.c:5188)()

fiona\_geometry.pyx in fiona._geometry.GeomBuilder.build (fiona/_geometry.c:3556)()

KeyError: 3002L

I have looked at the issue here and checked that 'id' 3002 exists, which it does as well as use Shapely to validate all of the geometries in the file.

I have another environment with the same versions of geopandas and a different version of fiona(1.5.1) that does work to load in the shp file, the one that doesn't work has fiona (1.7.0). Is it safe to say that this is a package issue with all this info? or am i missing something?


texts:
I'm looking for a cloud based database solution where I can store about 1 GB worth of tabular data and both myself and my client can access to perform SQL queries.  I'm looking for this database to have a simple Graphic User Interface that's intuitive to allow custom queries and to export results in a format such as CSV, DBF, etc. I'm okay with a monthly charge less than $50 per month.

I've looked at Amazon Web Services RDS, but I'm really confused how to "share" a MySQL database with my client. My understanding is that AWS RDS is a cloud based database solution where multiple users can access and manipulate data.  It's sounding more like a database has to be installed on your local machine and a connection is then made.  I'm not looking to download database software (ie MySQL), I just need the database and connections all to be in the cloud.


texts: ArcReader Publisher hung up in Exporting Raster Data

texts:
We recently upgraded to ArcMap 10.7 from 10.4. I never had a problem publishing my .mxd and creating a data package. Now, same .mxd gets stuck in the exporting raster data phase of creating a data package.

Any ideas that may have caused this?


texts: Writing a for loop in r to process multiple rasters with quantiles values

texts:
I am trying to get R to read a folder containing raster files, and then in a foor loop, extract quantile values from each of the raster to enventually process each raster with a substraction involving the quantile. The final step is to create a mean raster with the processed rasters but I'm struggling.

Here is what I've been writing so far:

path <- "C:/fichierR/rastermasque"

files <- list.files(path, pattern = "tif$", full.names = TRUE)

for(f in files){

decile <- quantile(f, probs = c(0.10))

r3 <- r2-decile

r4 <- calc(r3, fun = mean)

writeRaster(r4, "C:/fichierR/mean.tif", overwrite = TRUE)

}

The script when not in a for loop runs correctly but there are multiple errors in my attempt.


texts:
Is there a geopandas version of basemap's m.drawmapscale function?
If not, can I somehow utilize basemap's drawmapscale function within geopandas?


texts: Extracting individual images from gridded map

texts:
I am using QGIS (2.8.6-Wien) and I would like to do the following: Load a GeoTIFF, add a grid-layer (with cell size specified in pixels), select some cells and save each cell content as an individual image (e.g. as a .png file).

Since I am not very experienced with QGIS, I am not sure whether this is actually feasible with QGIS. What I managed to do so far is the griding. Below is a small example, where I used the Natural Earth Admin 0 - Countries dataset. I did the griding with "Vector" -> "Research Tools" -> "Vector Grid", Parameters X: 4.0, Y: 4.0 (output grid as polygons).

Now, I would like to select the cell indicated by "X" (Panama) and save the cell content as an image. I tried some of the QGIS geoalgorithms (Vector selection tools) but failed to reach my goal. Ultimately, I would like to select several cells and save each cell as an image, individually.

Could anyone give me a hint?


texts: How can I do client side linear referencing in an Esri Silverlight/WPF app?

texts:
I have an Esri Silverlight/WPF UI with a mapcontrol and a slider.  In the mapcontrol is a graphics layer with a single polyline on it, created using coordinates from a layer in a mapservice.  The shapefield of the featureclass used by the mapservice has PolylineM's.

I would like to find a way to move a point marker graphic along the polyline to reflect the M position as derived from the slider.  The min and max values of the slider are the MinM and MaxM values of the polyline.


texts: What is Esri's support plan for Python 2.7 beyond 2020 for ArcGIS Desktop?

texts:
According to PEP-373, Python 2.7 will retire on January 1, 2020 (Countdown Clock). That date was pushed back 5 years for those who couldn't migrate to Python 3.

According to this product life cycle page for ArcGIS Desktop, 10.6.1 will retire 4 years later on January 1, 2024. Sure, ArcGIS Pro supports Python 3, but there will be a gap of unsupported Python for those still using Desktop leading up to its retirement.

Esri's Deprecation Plan for 10.6 makes no mention of Python support beyond 2020.

I know that Python 2.7 is pretty stable, but I wonder, what are Esri's plans for support during this time?


texts: Model Only calculate value with zonal statistics as table

texts:
In a precedent post   Multiple Output for Zonal Statistics as Table I was trying to have as many table in the output as the rasters in the input. To overcome my errors, I decided to add a "Model Only Calculate Value" to the model. It looks like this:

I think it's not the proper way to connect the "calculate value" but i didn't know how to do it properly.
The error I have while running the model:

ERROR 000539: Error running expression: getFileName("C:\Users\fassi\Documents\Arcgis\NDVI_MODIS_per_decade")
Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 4, in getFileName
WindowsError: [Error 123] La syntaxe du nom de fichier, de répertoire ou de volume est incorrecte: 'C:\\Users\x0cassi\\Documents\\ArcGIS\\NDVI_MODIS_per_decade/*.*'

Failed to execute (Calculate Value).


texts: Reprojecting a shapefile into 'geographic coordinates' to load into iTree Canopy with QGIS

texts:
I have a polygon shapefile in CRS EPSG:27700 British National Grid that is a boundary of a study area that I would like to load into the iTree Canopy online tool.
https://canopy.itreetools.org/index.php

The tool uses Google Satellite as its basemapping so in QGIS I have saved the layer with CRS EPSG:3857 WGS84

But when I try to load it I receive the error 'Invalid projection format'

The pop-up help tells me that the shapefile must be 'reprojected to geographic coordinates'

I have looked at the help video which essentially describes what I have tried to do above using QGIS - ie. reproject the shapefile into WGS84.

I also tried based on another web search using the Vector > Geometry Tools > Export/Add Geometry Columns

Is there something obvious that I am not doing?


texts: Merging shapefile layers not working in QGIS?

texts:
How do i merge different vector layers in to one shapefile layer

In QGIS version 2.6


texts: How to group tools within a Python toolbox in ArcMap 10.2

texts:
I have created a simple Python toolbox at work for use in ArcMap 10.2. It is used by multiple members of staff so I'm hoping to make it as simple to use as possible. It currently only contains two scripts but I have others that I plan to add in the future. I am currently using the line below to load the two tools on startup, and whilst this is okay with only two scripts, it may become untidy as I add more.

self.tools = [opentemplatecsv, batchcsvtopdf]

Ideally it would be best to have one toolbox (i.e. "Dave's Tools") containing several sub-groups (e.g. "Mapping tools", "Modelling tools", "Hydrology tools"), each containing the appropriate scripts I've written.

Is there a way to place tools into sub-groups under the main toolbox, or is it best to simply add a few letters to the beginning of the scripts' names to alphabetically sort them?

i.e.

MAP - tool #1

MAP - tool #2

MOD - tool #3

HYD - tool #4

etc...


texts: How to get/pass token with ajax call to secured ArcGIS REST Service

texts:
I have a web map application that was written in HTML & JavaScript (using the ArcGIS JavaScript API). It desplays a few services that reside on an ArcGIS Server (10.2). It also has a custom legend that is created by requesting legend information directly from the ArcGIS REST Service using an ajax call.

It was recently decided that the application should be "locked down" and require a login in order to access it. We accomplished this by using the ArcGIS Server Manager to "lock" one of the services that the application consumes, and set a user/password to access it. The IdentityManager was also subsequently added to the application code. Now, when we navigate to the application, we are prompted for a username and password. After entering the correct credentials, access is provided and the "locked" layer successfully loads.

Unfortunately, the legend, which relies on an ajax call to the (now restricted) service endpoint, is not working. The ajax call returns error code 499 with the message "Token Required." I'm fairly new to ajax calls to begin with, and not at all familiar with service authentication and tokens. What is the correct way to pass the token with my ajax call in order to get the information I need to create my legend (and if I need to explicitely obtain the token, how do I do that?)


texts: Creating multiple points on Mapbox Isochrone map

texts:
Is it possible to create multiple points on an isochrone map in Mapbox?
This is the documentation for one of their maps with one point: https://docs.mapbox.com/help/tutorials/get-started-isochrone-api/
I want to do the exact same thing, just have multiple points with the buffers around them.


texts:
I am using the AHP add-on and ModelBuilder.

I am making the model in ArcGIS which deals with Multi-Criteria Decision Making Analysis (MCDM). For the part of assigning weights to criteria I am thinking to use AHP method. I have found AHP 2.0 add-on on codesharing.arcgis.com. It seems it works fine.
I have been thinking to include this in ModelBuilder. So when model comes to this part where AHP should be used, it will then open AHP 2.0 add-on. The user will work withing AHP as usuall and then once the criteria weights are caluclated it will be proceed further by taking these weights and caluclate the rest of MCDM.

How do I do this?

I anticpate that Python should be used here. I am familiar with this languge, but have not used it to establish the connection from one tool to another. In arcgis help website I have found the page - Integration external program within ModelBuilder. Here it explains the integration R program within ModelBuilder. I guess something similar has to be done in my case.

Do any of you have some tips, etc.. how this could be done?


texts: Creating vector line in the middle of two lines (grid) in QGIS 3.16

texts:
I have a grid. I want to create a new line in the middle of two lines. For example, calculate the distance between line 1 and 2 and create a new line (the red line in the image) exactly in the middle and of the same length as the grid lines, also between line 2 and 3 until the end of the grid.
The attribute table contains the number of rows and their type.


texts: Select all images in image collection Google Earth Engine in R Rgee to mosaic

texts:
Is there a way to select all images in an image collection in Google Earth Engine or using Rgee?
I want to select all the images in an image collection after filtering the collection and then mosaic the images.
I know I can select individal images this way:
s2_sr16<-ee$ImageCollection('COPERNICUS/S2_SR')$filterBounds(SHP)$filterDate('2016-01-01', '2016-12-31')
listOfImages <- s2_sr16$toList(s2_sr16$size())
selected_images <- ee$ImageCollection$fromImages(list(ee$Image(listOfImages$get(1)), ee$Image(listOfImages$get(2)),
                                                      ee$Image(listOfImages$get(3)), ee$Image(listOfImages$get(4)),
                                                      ee$Image(listOfImages$get(5)), ee$Image(listOfImages$get(6))))


texts: gdalwarp leaves big hole at 180° longitude when warping from 4326 to custom AEqD

texts:
I'm trying to reproject the Natural Earth I 1:50m raster into an Azimuthal Equidistant projection centered on my hometown. As you may know, the source raster covers from -180° to 180° latitude, and -90° to 90° Here's the command I'm using:

gdalwarp -multi -dstnodata 0 -ot Byte -t_srs "+proj=aeqd +lat_0=48.037 +lon_0=-122.408 +ellps=sphere +units=m +no_defs" -r cubic -of GTiff -ts 10800 10800 NE1_50M_SR_W.tif natural_earth_10800.tif

This works pretty OK, except that it leaves a big, ugly gap west of 180°. Also, it cuts off some of each side of the output raster.

I've experimented with setting a radius instead of the +ellps parameter, with no change. Using WGS84 gives only a small slice of the Earth as output.

Here's what I got when I tried +R_A instead of +ellps=sphere:

The result with WGS84 was much the same. Also tried reprojecting to Spherical Mercator as an intermediate projection before reprojecting to my custom AEqD, same artifacts appear.

So, how do I fill in the gap? The pixels are definitely there in the input raster, they're just not making it into the output.


texts: Buffer option not available to select in Editor menu of ArcMap?

texts:
I need to use the Buffer in the Editor menu of ArcMap, but it is unavailable.

What is the problem?


texts: How to mask a raster layer to only include specific values in QGIS?

texts: I'm attempting to reduce a raster layer, wherein each point represents night-time luminosity/intensity of city lights, to include only those pixels with values in the top 10% of values. Effectively, I'm hoping to generate a new raster layer that only includes the brightest points. I'd imagine I could do so via some sort of 90% quintile masking specification but am at a loss as to how this might be accomplished in QGIS. Doing so via the expression table for a vector layer makes sense to me, but how might this be done in the case of a raster?


texts: Creating python dictionary that maps each domain to multiple coded values in File Geodatabase?

texts:
I have a File GDB that has domains and coded values. I would like to extract those values from the GDB and build a dictionary. For example:

I have hundreds of domains and with each domain has a set of coded values associated with it.

*Domain* , *Coded Value*
['tree', 'pine']
['tree', 'magnolia']
['tree', 'fir']
['tree', 'grey pine']
['tree', 'oak']
['soil', 'clay']
['soil', 'loam']
['soil', 'sandy']
['brush', 'manz']
['brush', 'short']
['water', 'wetland']
['water', 'lake']
etc.....

I want this:

{'tree': ['pine', 'magnolia', 'fir', 'grey pine']}
{'soil': ['clay', 'loam', 'sandy']}
{'brush': ['manz', 'short']}
{'water': ['wetland', 'lake']}

How do i accomplish this by using arcpy.da.ListDomains? This is what I have so far:

import arcpy

doms = arcpy.da.ListDomains(gdb)

    for dom in doms:

    if dom.domainType == 'CodedValue':
        codedvalues = dom.codedValues
        for code1 in codedvalues:

             list1 = []

             c_data = "{},{}".format(dom.name, code1)
             domain2 = c_data.split(",")[0]
             code2 = c_data.split(",")[1]

             list1 = [domain2, code2]
             print list1  #this prints out the first code block shown above


texts: Importing legacy shapefiles and attribute information to geodatabase

texts:
I am trying to take older sewer/water system shapefile data with limited, but important attribute information and importing it in to a feature class in a geodatabase, while keeping the attribute information. These older files do not have similar field names as the geodatabase in which to import the attribute information.

Any advice in which to complete this task is greatly appreciated. The solution does not neccesarily have to be solved programmatically.


texts: Can the correct TIGER Line Id be returned for an address from the PostGIS Tiger Geocoder?

texts:
After reviewing the PostGIS Tiger Geocoder documentation, I don't see a method of returning the TIGER Line Id for the address.  I would assume this is information is available internally as part of the the street range interpolation processing.

http://postgis.net/docs/Extras.html


texts: How large must be my raster for rasterizing a shapefile in R in a scale 1:100.000?

texts:
I´m trying to rasterize a shapefile of land cover for the colombian coast. The original shapefile is in a scale of 1:100000. I´m using this code in R language.

costa <- readOGR(file.choose()) #Here, I put the original shape
r <- raster(ncol=900000, nrow=900000) #This is the line in which I have problems. How big should be my raster for not losing too much information?
extent(r) <- extent(costa)
rcosta <- rasterize(costa, r)


texts: Determine if a raster layer exists in a file geodatabase

texts:
I'm trying to determine if a layer file already exists in a given file geodatabase using arcpy.Exists.

I wrote a code that works without problems for Feature layers, however it doesn't work for Raster layers.

To figure out where the problem is, in my code I create a new empty gdb and check if a layer baseName exists in the new gdb. For Feature layers it returns False, as one would expect, given that the gdb is empty. However it returns True for the raster layers. The code is exactly the same except one looks at Feature layers and the other at Raster layers. Where is my mistake?

Here is the code for Feature layers, returning False for all Feature layers:

import arcpy

 #sets the processing target to the current aprx

aprx = arcpy.mp.ArcGISProject(r"C:/Test.aprx")

# Execute CreateFileGDB
arcpy.CreateFileGDB_management("C:/empty.gdb")

m = aprx.listMaps("Map")[0]
#loop to gather source paths of each layer in current aprx and export each to the out directory defined above

arcpy.env.workspace = r"C:/empty.gdb"
for lyr in m.listLayers():
#checks that each layer has a datasource attribute
    if lyr.supports("DATASOURCE") and lyr.isFeatureLayer is True:

       #grabs the layer's datasource
        ds = lyr.dataSource
        # Describe the feature class
        desc = arcpy.Describe(lyr)
        #grabs the layer's name from the table of contents
        nm = desc.baseName
        print(arcpy.Exists(nm))

Here is the code for Raster layers, returning True for all Raster layers:

import arcpy

 #sets the processing target to the current aprx

aprx = arcpy.mp.ArcGISProject(r"C:/Test.aprx")

# Execute CreateFileGDB
arcpy.CreateFileGDB_management("C:/empty.gdb")

m = aprx.listMaps("Map")[0]
#loop to gather source paths of each layer in current aprx and export each to the out directory defined above

arcpy.env.workspace = r"C:/empty.gdb"
for lyr in m.listLayers():
#checks that each layer has a datasource attribute
    if lyr.supports("DATASOURCE") and lyr.isRasterLayer is True:

       #grabs the layer's datasource
        ds = lyr.dataSource
        # Describe the feature class
        desc = arcpy.Describe(lyr)
        #grabs the layer's name from the table of contents
        nm = desc.baseName
        print(arcpy.Exists(nm))


texts: "CountIf" or alternative function in ArcMap field calculator?

texts:
Within ArcMap 10.3 field calculator, I am looking to count the number of features in a column that match the current value, and output that count. Ideally, the output would look like this:

Number       Count
2000              2 (because there are only 2 counts of 2000 in the number column)
2000              2 (because there are only 2 counts of 2000 in the number column)
2001              1 (because there is only 1 count of 2001 in the number column)
4333              1 (because there is only 1 count of 4333 in the number column)

I accomplished this in excel using the Visual Basic Countif() function, but I just don't know how to call a full array of values from the Number column.


texts: Save separate attribute values using Relation reference in QGIS?

texts:
I have a list with different possible (allowed) combinations of different attributes (values) in a csv-file.
For example vegetation.csv

Number;Letter;Code;Bush;Soil;Tree
1;a;1a;k;x;f
1;a;1a;k;y;
1;a;1a;l;v;g
1;b;1b;;p;h
1;b;1b;m;q;h
1;b;1b;m;r;
2;a;2a;s;;f
2;b;2b;;x;
2;c;2c;;l;h
2;d;2d;k;;h
2;d;2d;k;v;g
2;e;2e;w;s;f
3;a;3a;b;;g
4;a;4a;b;v;h

I manage to get "depending" combo boxes by using project relations and the Relation reference widget with filters and chain filters: If I choose Number = 1, I may choose Letter a or b. If I choose b, I may set Bush to empty or m. If m, I can set Soil to q or r. If r, Tree is set to empty.

I need to save 4 attributes: Code, Bush, Soil and Tree.
so far I only succeed in saving Code (since Code is the referred field in the Project relation setup).
Is it somehow possible to save all the 4 attributes "in one operation"?

I may create more project relations, but then I will have to repeat the combo boxes 4 times (one time for each attribute I want to store)? I want the user to only once chose the combination of Number, Letter, Bush, Soil and Tree, and that these choses are saved into the 4 attributes Code, Bush, Soil and Tree.

Is this achievable?


texts: How to identify feature from multiple layers in QGIS?

texts:
I have recently migrated from ArcGIS to QGIS.

Is there any way to identify a feature having multiple layers? I am currently working with about 400 shape files and need to know which is the one I'm currently working on. When I click the polygon I want to identify with the identify tool nothing happens unless I have the layer selected in the layer panel. The same happens with the selection tool.

How can I fix this? In ArcGis I could identify any feature without having it selected on the layer list and I've looked through the settings of QGis and found nothing.


texts: How to divide a polygon into many polygons based on overlapping layer

texts:
I am new to GIS and also QGIS tool.

There are two layers imported into the qgis project.
Layer-1 (green) has attributes ID, area, code
Layer-2 (brown) has attributes station, shp_area, name.

I want to divide the polygons by overlapping and then add a new field value with the divided polygon's respective area. Other fields will remain the same of Layer-1
For e.g. Along the border of one polygon (brown border), a larger polygon (green) was divided into two polygons. The total area before division was 125sqm. After division (shown in red), the areas respectively became 70sqm and 55sqm.

I want to do this in the entire shapefile and store the new output (divided green areas) as Red-layer.


texts: OpenLayers - createBox() relative to map view rotation

texts:
I would like to make ol.interaction.Draw.createBox() method to draw the box relatively to map view rotation. By default it is aligned with the coordinate system axes.
Current behavior:

Desired behavior:

I have tried making custom geometryfunction for ol.interaction.Draw to calculate the coordinates of the box with map view rotation taken into consideration but cannot get it to work.
It draws weird shapes like below:

My geometryfunction is as follow:
var geometryFunction = function(coordinates, geometry) {

    if (!geometry) {
    geometry = new ol.geom.Polygon(null);
    }

    var currentViewRotation = map.getView().getRotation();

    var firstCorner = coordinates[0];
    var thirdCorner = coordinates[1];

    var dy = thirdCorner[1] - firstCorner[1];

    var secondCornerX = firstCorner[0] - dy * Math.sin(currentViewRotation);
    var secondCornerY = firstCorner[1] + dy * Math.cos(currentViewRotation);
    var forthCornerX = thirdCorner[0] + dy * Math.sin(currentViewRotation);
    var forthCornerY = thirdCorner[1] - dy * Math.cos(currentViewRotation);

    var secondCorner = [secondCornerX, secondCornerY];
    var forthCorner = [forthCornerX, forthCornerY];
    var fifthCorner = coordinates[0];

    var newCoordinates = [];
    newCoordinates.push(firstCorner);
    newCoordinates.push(secondCorner);
    newCoordinates.push(thirdCorner);
    newCoordinates.push(forthCorner);
    newCoordinates.push(fifthCorner);
    geometry.setCoordinates([newCoordinates]);

    return geometry;

};

I have tried playing with the sin/cos and +/- in calculating the coordinates but nothing worked.
Does anyone have this figured out?


texts:
I have a feature class with about 250,000 contiguous polygons.  ArcEditor turns up around 9,000 overlap errors when I validate the topology.  I want to merge each of these errors into one of the polygons they source from - they're all very tiny overlaps, so it really doesn't matter which one.  All 9,000 errors can be fixed in exactly the same way.

In the Error Inspector, the Merge option is only available when I select one error at a time - selecting a bunch of them greys it out.  This means I need to do the same repetitive task for every single error.  I would love to be able to avoid doing these one at a time, might save me from carpal tunnel and wasting the better part of a day.

Any suggestions?


texts: How to add/upload own vector file (Shape, KML, GeoJSON) to Openlayers and show it on a map

texts:
I started working with Openlayers and web mapping. Therefore I use a Cartaro (a web mapping platform) based on Drupal. Because I'm not a programmer it is (at the moment) much easier for me to administrate OpenLayers with the graphical interface in Drupal/Cartaro. Now I would like to add the following feature and I think it is only possible with some programming effort in OpenLayers:

I would like to give users the chance to upload their own vector file (mostly Shapefile, KML or GeoJSON) and show it temporarily on the map.

Does anyone of you has an idea or an example how to do this?

Thanks in advance!


texts: Geoserver "Could not find a style for layer" error

texts:
I have a bunch of layers stored in Geoserver. However, when I try to preview certain layers, a WMS exception gets thrown.

<ServiceExceptionReport version="1.1.1" >
    <ServiceException code="NoDefaultStyle">
        Could not find a style for layer wyoming:GoatHuntAreas2014, either none was specified or no default style is available for it
    </ServiceException>
</ServiceExceptionReport>

When trying to manage the layer, I click on the "Publishing" tab and I see the following error message (truncated).

Any ideas? How can I go back and set a "default style?"

org.apache.wicket.WicketRuntimeException: Exception in rendering component: [MarkupContainer [Component id = defaultStyle]]
at org.apache.wicket.Component.renderComponent(Component.java:2725)
at org.apache.wicket.MarkupContainer.onRender(MarkupContainer.java:1538)
at org.apache.wicket.Component.render(Component.java:2517)
at org.apache.wicket.MarkupContainer.renderNext(MarkupContainer.java:1440)
at org.apache.wicket.MarkupContainer.renderComponentTagBody(MarkupContainer.java:1603)
at org.apache.wicket.MarkupContainer.onComponentTagBody(MarkupContainer.java:1527)
at org.apache.wicket.Component.renderComponent(Component.java:2686)
at org.apache.wicket.MarkupContainer.onRender(MarkupContainer.java:1538)
at org.apache.wicket.Component.render(Component.java:2517)
at org.apache.wicket.MarkupContainer.renderNext(MarkupContainer.java:1440)
at org.apache.wicket.MarkupContainer.renderComponentTagBody(MarkupContainer.java:1603)


texts: in cursor.updateRow(row) RuntimeError: Object invalid or no longer set

texts:
I'm running an arcpy code in IDLE that update a field in shapefile. After updating it get the following error:

  Traceback (most recent call last): File "C:/Users/10/Desktop/test.py", line 21, in cursor.updateRow(row_st) RuntimeError: Object invalid or no longer set.

And, here is the code:

from arcpy.sa import *
import arcpy
from arcpy import env
env.workspace = r"D:\Phd-Files\dataset.mdb"
env.overwriteOutput = True

fields_Street=["OBJECTID","SHAPE@","t_light_de", "st_dense", "bikeLn_den", "Buff_area", "len"]

cursor= arcpy.UpdateCursor("Boston_Street_Buff_10_Joint", fields_Street)
for row_st in cursor:
    geom = row_st.getValue('Shape')
    # calculate street density
    _area_buff=0
    arcpy.Buffer_analysis(geom, "geom_buff", "10 Meters")
    cursor_buff= arcpy.SearchCursor("geom_buff", ["ORIG_FID", "Shape_Area"])
    for row_buff in cursor_buff:
        _area_buff=_area_buff+ row_buff.getValue('Shape_Area')
        row_st.setValue ('st_dense' , row_st.getValue('Buff_area')/_area_buff)
    del cursor_buff, row_buff
    print ('street_ObjectID={0}, bike_lane density= {1}'.format(row_st.getValue('OBJECTID'), row_st.getValue('st_dense')))
    cursor.updateRow(row_st)
    arcpy.Delete_management("geom_buff")


texts: Invalid multipolygon of valid individual polygons

texts:
Could anyone please explain this to me? The polygons all have valid geometry but the resulting multipolygon somehow doesn't. The following script:

from shapely.geometry import Polygon, MultiPolygon

polygon1 = Polygon([(0,0), (1,0), (1,1), (0,1), (0,0)])
polygon2 = Polygon([(1,1), (2,1), (2,2), (1,2), (1,1)])
polygon3 = Polygon([(1,0), (2,0), (2,1), (1,1), (1,0)])
polygon_list = [polygon1,polygon2, polygon3]

mpolygon = MultiPolygon(polygon_list)

for item in polygon_list:
    print(item.is_valid)
for item in mpolygon:
    print(item.is_valid)

print(mpolygon.is_valid)

returns the following:

True
True
True
True
True
True
False

I believe this is causing me problems as I need to intersect a multilinestring with a multipolygon, and the intersection method requires for all the polygon geometries to be valid.

Why are the individual polygons valid while the full multipolygon is not?


texts: Is it possible to create map packages for ArcMap version 10 using 10.1?

texts:
Is it possible to create a map package in an older version? Currently I am using 10.1 and the client uses 10.0. Can you create a map package for different versions?

I see you can create layer packages using Python for different versions but you cannot specify the version when creating a map package. As I don't have 10.0 I can't test whatever I create these.

When you unpack the file using ArcCatalog does it create it v.10. and 10.1?

Can you create the original file in 10.0?


texts:
I need to build a website which needs to perform the following tasks:

Display a web map with topographical details
On the map show a set of selected streets
When a street is selected show a set of points of interest on that street

The constraints are:

Ruby on Rails as a programming language

Everything in the stack must be open source
I've given an interested look at OpenStreetMap will it work for the needed tasks?

EDIT
Actually there is a product using Google Earth in this company but they want to change it due to its commercial cost.


texts: How to do Aerial Triangulation from raw aerials?

texts:
I have received 1390 raw aerial images. For the purpose of understanding I have imported 4 images and here you will see that roads are not lining up.

When I spoke to the creator of these aerials, he mentioned that inorder to match those  I need to do aerial triangulation. They have also sent me Exterior Orientation Computation Utility file having fields:

ID, # EVENT, TIME (s), EASTING, NORTHING, ORTHOMETRIC HEIGHT, OMEGA, PHI, KAPPA, LAT, LONG.

How do I perform the aerial triangulation in an open source tool or any other paid version program?


texts: Deleting feature dataset containing network dataset from ArcSDE database?

texts:
I have trouble with deleting unwanted feature datasets from ArcSDE.

Datasets contain network datasets which cannot be deleted either and I think they cause the problem. Options to delete or remove are simply greyed out. I haven't found anything useful in the web so far about this, any suggestions?


texts:
I am trying to apply the style from .qml to a group of layers in QGIS. I have tried to incorporate LoadNamedStyle within jkall's code here:
https://github.com/jkall/qgis-processing-scripts/blob/master/set_transparency_for_group.py
I end up with this:

from qgis.core import *
def applystyle_group(name):
    root = QgsProject.instance().layerTreeRoot()
    point = root.findGroup(name) #Find Group AP
    for child in point.children():
        if isinstance(child, QgsLayerTreeLayer):
            if child.layer().type()==0:
                child.loadNamedStyle('*.qml')
applystyle_group("AK")

It gives me an error:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "C:/OSGeo4W64/bin/applystyle_group.py", line 9, in <module>
     applystyle_group("AK")
  File "C:/OSGeo4W64/bin/applystyle_group.py", line 8, in applystyle_group
    child.loadNamedStyle('*.qml')
AttributeError: 'QgsLayerTreeLayer' object has no attribute 'loadNamedStyle'

Does QgsLayerTreeLayer have something similar to loadNamedStyle?


texts:
I've got a polygon dataset with a few different variables.

I want to be able to see if these variables are spatially auto-correlated and then also run a Geographically Weighted Regression for the variables.

I feel like I did this in ArcGIS a long time ago and it was fine, but when I look at QGIS' tools for these things, I can only seem to input point datasets?

Is there something I am missing?


texts: Adding scaled cost line barriers to Network Analyst > Closet Facility

texts:
I am trying to use an existing line shapefile to create a preferred route using scaled cost line barrier in network analyst closet facility tool. I thought by using the existing line file as a scaled cost line barrier with Attr_cost set to less than 1 it would be more likely for the closet facility to route using that path over the network dataset route from streetmaps routing network dataset.
I may be approaching this incorrectly but I'm trying to use our existing network as the main routing tool and only use the streetmap routing secondarily. In a sense existing versus new build network.


texts:
I have drawn a series of straight lines on a map joining different locations. The lines have a single vertex at each end.  I am using the Ordnance Survey of Great Britain 1936 projection.

I now want to know the compass bearings I will have to walk to move from one location to the next.

Is there any easy way to get QGIS to calculate/display the azimuth of the lines I have drawn?
Better still is there any way to get QIS to actively display the azimuth of the lines as I am drawing them?

I am aware there is an azimuth function but the syntax appears to need the input of the coordinates of the vertices on either end of the line which would be very time consuming.


texts: How to improve Triangulated Irregular Network (TIN) display speed?

texts:
I am preparing Triangulated Irregular Network for large area (5000 Sq Km) in ArcGIS, The TIN models takes too much time to open in ArcMap, and if I move or change anything again its starts to rendering, It takes too much time to do small changes, how to stop the rendering, and how to keep it in Builded TIN in ArcMap.

Number of Nodes in TIN is 3633846
Number of Triangles is 7267645

it takes three minutes to open fully.


texts: RTK Coordinates from Phone to WGS84 Dec Deg using ITRF00 Transformation and proj4.js

texts:
I am attempting to simulate an RTK GPS unit transmitting a CORS corrected GPS location to a phone. The coordinates are in NAD83 (from CORS) but in the phone they are seen in decimal degrees. I need to get the coordinate from NAD83 to WGS84 in order to properly place the location on the map (which uses the WGS84 datum). I want to use the ITRF00 transformation, with the definitions coming from the latest ArcGIS Desktop 10.7 and adjusted with input from here.

The problem I am having is when I do the transformation the output is not in decimal degrees, and if I convert back to wgs84 with no transformation, the result is not transformed.

var nad83 = "+proj=longLat +ellps=GRS80 +datum=NAD83 +towgs84=-0.9956,1.9013,-0.5215,0.025915,0.009246,0.011599,-0.00062 +units=degrees +no_defs"
var nad83NoTransform = "+proj=longLat +ellps=GRS80 +datum=NAD83 +towgs84=0,0,0 +units=degrees +no_defs"

Here are my efforts -

https://codepen.io/getbounds/full/BaaEVKy

If there is another library that can do this let me know. This needs to be done in the browser.


texts:
I am wondering how to find the parent map tiles of a particular tile or the tiles contained within a particular map tile. I know tiles usually have the structure of somedir/{z}/{x}/{y}.png/jpg and I know their is a formula for projecting WGS84-coordinates to tile-names:

http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames

E.g. If you change a geometry on a zoom level 18 map tile, do you then have to calculate the tilename on zoomlevels 17 - 1 using the formula to find out which tiles are effected or is there a more efficient way of finding the parents or children of a certain tile?

I have also found this link with notes on how to update tiles, but there seems to be no information concerning my question:

http://wiki.openstreetmap.org/wiki/Tirex/Tile_Update_Strategies


texts: QGIS 3 system date format is not assumed

texts: The system date format is 'dd/MM/yyyy', and QGIS recognizes this on the Options / General tab, however it does not assume this format anywhere else in itself. Any loaded feature that has a date or datetime attribute has the format 'yyyy-MM-dd'. How to set QGIS 3 default date format so that it actually uses the system format?


texts: Error code '500' : 'Error performing submit job operation

texts:
I have an error in line which calling SubmitAndPollStatusAsync(parameter) Function

    public async void DisplayGrid()
    {
        try
        {
            Uri GridUri = new Uri("http://rsarcgis.eastus.cloudapp.azure.com/arcgis/rest/services/Geoprocessing/GenerateGrid/GPServer/GenerateMapGrid");

            _gpTask = new Geoprocessor(GridUri);

            var parameter = new GPInputParameter();
    parameter.GPParameters.Add(new GPLinearUnit("Width", LinearUnits.Miles, Width));
            parameter.GPParameters.Add(new GPLinearUnit("Height", LinearUnits.Miles, Height));
            parameter.GPParameters.Add(new GPLong("Rows_Count", RowsCount));
            parameter.GPParameters.Add(new GPLong("Columns_Count", ColumnsCount));
            parameter.GPParameters.Add(new GPString("Origin_Coordinate", CoordinateX+" "+CoordinateY.ToString()));

            parameter.OutSpatialReference = Globals.AppMainMapView.SpatialReference;
            parameter.ProcessSpatialReference = Globals.AppMainMapView.SpatialReference;
            var result = await SubmitAndPollStatusAsync(parameter);
            if (result.JobStatus == GPJobStatus.Succeeded)
            {

                var resultData = await _gpTask.GetResultDataAsync(result.JobID, "Clipped_Counties");
                if (resultData is GPFeatureRecordSetLayer)
                {
                    GPFeatureRecordSetLayer gpLayer = resultData as GPFeatureRecordSetLayer;
                    if (gpLayer.FeatureSet.Features.Count == 0)
                    {
                        var resultImageLayer = await _gpTask.GetResultImageLayerAsync(result.JobID, "Clipped_Counties");

                    }
                }
            }
        }
        catch (Esri.ArcGISRuntime.Http.ArcGISWebException exp)
        {
            Messenger.Default.Send(exp);
        }
    }

Code of SubmitAndPollStatusAsync( ) Function is :

      private async Task<GPJobInfo> SubmitAndPollStatusAsync(GPInputParameter parameter)
    {
        // Submit gp service job
        var result = await _gpTask.SubmitJobAsync(parameter);

        // Poll for the results async
        while (result.JobStatus != GPJobStatus.Cancelled && result.JobStatus != GPJobStatus.Deleted
            && result.JobStatus != GPJobStatus.Succeeded && result.JobStatus != GPJobStatus.TimedOut)
        {
            result = await _gpTask.CheckJobStatusAsync(result.JobID);

           // txtStatus.Text = string.Join(Environment.NewLine, result.Messages.Select(x => x.Description));

            await Task.Delay(2000);
        }

        return result;
    }


texts: Obtain High Resolution Graphics suitable to Print

texts: I have created a map with a single layer that needs to be exported for printing (30 dpi quality and quite large). I researched this extensively and see that the ability to export an image has been deprecated in the new Carto Builder UI. How should I proceed with getting this exported?


texts: Getting aspect and reclassifying raster files using Python in a loop

texts:
By using Python in PythonWin 2.7.3, I tried to obtain aspect files of a list of raster files in a folder(called database_Dem20) and save it to another folder, database_Dem20/testOutput.
In the same loop, after getting the aspect file, i want to reclassify this aspect file by using reclassify toolbox.

The first part, getting aspect, is doing fine. However, the reclassify part has error, ERROR 000622: Failed to execute (Reclassify). Parameters are not valid.

I have read the syntax from Esri website for many times and tried to type in different ways but still couldn't figure out where is wrong.

Can anyone help me, please?

Here is my Python scripts:

# Import arcpy module
import arcpy
from arcpy import env
import os
from arcpy.sa import *

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True

# Set Geoprocessing environments
env.workspace = "C:\\database_Dem20"
targetWorkspace = "C:\\database_Dem20\\testOutput"
print "==============================NEW START================================"

# Process: Aspect
# Get and print a list of TIFs from the workspace
rasters = arcpy.ListRasters("*", "TIF")

for raster in rasters:
    inRaster = raster
    outAspect = Aspect(inRaster)
    outAspect.save(targetWorkspace + "\\AA_" + inRaster)

    print("------------------------Reclassify---------------------------------------")
    # Process: Reclassify
    inAspect = str(outAspect)
    remap = RemapRange([[-1, -9.9999999999999995e-007, 1],[-9.9999999999999995e-007, 22.5, 2],[22.5, 67.5, 3],[67.5, 112.5, 4],[112.5, 157.5, 5],[157.5, 202.5, 6],[202.5, 247.5, 7],[247.5, 292.5, 8],[292.5, 337.5, 9],[337.5, 360, 10]])
    outReclassify = Reclassify(inAspect, "Value", remap ,"NODATA")
    outReclassify.save(targetWorkspace + "\\BB_" + inAspect)
    print("-------------------------Reclassify--------------------------------------")
del raster
del rasters


texts: Extract filename from path in field calculator

texts:
I try to extract the filename from a path in my attribute table in QGIS.

regex_substr should be the way to go I assume, but testing all kinds of example from the web did not work out.

I have images which I want to link, but the original pathes are not valid anymore, cause all images from different people went into one single folder now.

Example paths:

./20170608-attachments/p__20170608_122128.jpg
./20170603-attachments/p__20170603_095124.jpg

here I would just like to populate a new field with p__20170608_122128.jpg etc. only.

A solution, that would give me the rest of the string after the last / would be the most versatile one in my case, I assume.

Any ideas?


texts:
I have multiple layers with 'many features' within. One layer is a larger boundary layer which contains these 'many features'. The example below shows my situation. The red outline polygon is the boundary layer and the various coloured layers within are the 'many features' which have an all important AREA (m2) field in each. For line features, this is LENGTH (m).

Ideally I want to spatially join all the 'many features' to each Boundary layer using 'HAVE_THEIR_CENTER_IN'.
However, this spatial join method looks at the centre of the Boundary polygon, not the centers of the multiple features. ie. It will join to the boudary layer, any layer which sits within the center of the boundary layer...
I want it to join any feature that has their center within the boundary layer, TO the boundary layer.

how should I best work around this problem?
It is the wrong way round for my task.


texts:
I am wondering if there is an optimum raster tile size when using -t in raster2pgsql. I am setting -t auto and it creates thousands of tiles and takes hours to load.

I am working with very large rasters with 1m resolution.

I know that there is a limit of ~55,000 rows X ~55,000 columns. Should this be my tile size? Does auto produce a raster that will be processed the fastest in subsequent analysis using, e.g., map algebra?


texts: How to convert .HDF to .Tiff and keep the number of the output files exactly as the input?

texts:
I have 332 .HDF files that I have processed in R to convert them to .Tiff. It has worked fine, however, the number of output files are not the same as the number of the input files! Input is 332, output is 328. Does anyone knows why this is happening and how it can be resolved, please?

This is the code i have used for converting:

for (i in 1:length(hdfs)){

setwd("F:\\MODIS\\MOD13Q1.005\\h21v08\\")

reprojectHDF(hdfName       = hdfs[i],
           filename      = paste(d_temp,"\\",out_names[i],sep = ""),
           MRTpath       =  MRTpath,
           resample_type = "NEAREST_NEIGHBOR",
           proj_type     = "GEO",
           bands_subset  = "0 0 1 0 0 0 0 0 0 0 0 0 0",
           proj_params   = "0 0 0 0 0 0 0 0 0 0 0 0",
           datum         = "WGS84",
           pixel_size    = 0.004514)}


texts: How can I determine elevation levels for a region in Java?

texts:
I need to find out the elevation levels for a certain region (center is located at 132.159633, 43.350116).

For this purpose I went to the Elevation finder and downloaded there the relevant part of GTOPO30 data set.

Now I have two files with elevation data:

gt30e100n90.tif
Contents of the gt30e100n90_dem.zip archive, incl. gt30e100n90.dem.

How can I extract the elevation values (input: lon/lat coordinate, output: elevation at this location according to those data) in Java?

Note that I will only retrieve elevation data for buildings, hence I don't need to take into account that a particular coordinate may correspond to water.

I saw this project, there seems to be no source code there.


texts:
Zurich 3.81 Imported photo with metafile into project. Set display label to show Attributes Lat and Lon coordinates. Set the decimal places to 8. All worked fine except some of the Lat/Lon fields when imported resulted in recurring digits which the set decimal places did not recognise. Despite saving layer/project/rebooting nothing worked apart from manually cutting these aberrant cells back to 8 places.
Is this a bug or have I failed to apply something correctly?
import


texts: How to create custom color ramp using "r.colors" GRASS GIS function?

texts:
I am trying to assign a raster/grid a custom color palette in GRASS using the "r.colors" function.

r.colors map=Consensus_BRT_RCP85 colors=rules

This allows me to input my rules for my color palette:

0 238:238:238
0.06 0:0:255
0.33 0:255:255
0.67 255:255:0
1.0 255:0:0
end

This function operates, but it doesn't create the color palette I would like.

I'd like it to be like my custom GMT color palette:

0.000000 238 238 238 0.060000 238 238 238
0.060000   0   0 255 0.333333   0   0 255
0.333333   0 255 255 0.666667 255 255   0
0.666667 255 255   0 1.000000 255   0   0

My GMT palette looks like this (I want my probabilistic map to have this color ramp):

I've seen that there is a function called "r.cpt2grass" but this function isn't found on my GRASS installation (not sure what to do about GRASS GIS "add-ons").


texts: Removing Clouds From (Landsat8) NDVI Imagery with ArcMap?

texts:
I recently was able to generate NDVI images from Landsat 8 using Arcmap. My end goal is to run some zonal statistics over specific neighborhood polygons, however with the cloud values (ranges of 0.002< X <0.025) still in my images, this will skew my results. I've thus far tried to figure out how to remove this using raster calculator but I'm not certain I'm actually doing what I think I am.

Using ArcGIS Desktop 10.5.1, how could I remove the above data range without altering the rest of my image?


texts: Turning Polyline with Multiple Segments into One Single Line in QGIS?

texts:
Using QGIS Valmiera (MapInfo and other software is also available)

Is there a way to convert a polyline with multiple segments into a single line? I.e. by taking the first point and the last point of the polyline and simply connecting the two.

Furthermore, is this possible to be done under a condition that say the polyline is straight enough to be considered a single line. For example if polyline is within a certain distance of the line connecting the start and the end point.

The context is I'm importing into a hydraulics software package with a pipe limit and polylines are disaggregated in this software and as a result count for thousands of pipes rather than hundreds despite being largely straight.


texts: Change symbols outline color without losing fill pattern

texts:
I have defined a layer symbology like shown in the screenshot below (categories are applicable laws however).

The layer symbology is defined like this i.e. (no pattern fill):

For these three categories it would be no big deal to change the outline color separately, but imagine some tens of categories this would be linear time consuming.

So is there a way to change some aspect of the symbology (like the outline color) at once without losing the settings for the fill style?


texts: Moving marker along path/polyline using Leaflet API?

texts:
Trying to migrate from Googlemaps to Leaflet, some errors return in following functions.

Number.prototype.toRad = function() {
  return this * Math.PI / 180;
}
//
Number.prototype.toDeg = function() {
  return this * 180 / Math.PI;
}
// polyline length calc
L.Polyline.prototype.m_len = function () {
  var distance = 0,  coords = null, coordsArray = this._latlngs;
  for (i = 0; i < coordsArray.length - 1; i++) {
    coords = coordsArray[i];
    distance += coords.distanceTo(coordsArray[i + 1]);
  }
  return distance;
}
//
function moveAlongPath(point, distance, index) {
  index = index || 0; // Set index to 0 by default.
  if (index < routePoints.length) {
    var poly = new L.Polyline([routePoints[index], routePoints[index + 1]]);
    // alert(JSON.stringify(poly, null, 4));
    var distanceToNextPoint = poly.m_len();
    if (distance <= distanceToNextPoint) {
      return routePoints[index].moveTowards(routePoints[index + 1], distance);
      // error: routePoints[index].moveTowards is not a function
    } else {
      var d_distance = distance - distanceToNextPoint;
      return moveAlongPath(routePoints, d_distance, (index + 1));
    }
  } else {
    return null;
  }
} // end moveAlongPath
/////////////////////////////////
L.LatLng.prototype.moveTowards = function(point, distance) {
  var lat1 = this.lat.toRad();
  var lon1 = this.lng.toRad();
  var lat2 = point.lat().toRad();
  var lon2 = point.lng().toRad();
  var dLon = (point.lng() - this.lng()).toRad();
  // Find the bearing from this point to the next.
  var brng = Math.atan2(Math.sin(dLon) * Math.cos(lat2),
                        Math.cos(lat1) * Math.sin(lat2) -
                        Math.sin(lat1) * Math.cos(lat2) *
                        Math.cos(dLon));
  var angDist = distance / 6371000 ; // 66371000 Earth's radius.
  // Calculate the destination point, given the source and bearing.
  lat2 = Math.asin(Math.sin(lat1) * Math.cos(angDist) +
                   Math.cos(lat1) * Math.sin(angDist) *
                   Math.cos(brng));
  lon2 = lon1 + Math.atan2(Math.sin(brng) * Math.sin(angDist) *
                           Math.cos(lat1),
                           Math.cos(angDist) - Math.sin(lat1) *
                           Math.sin(lat2));
  if (isNaN(lat2) || isNaN(lon2)) {return null; }
  return L.LatLng(lat2.toDeg(), lon2.toDeg());
} // end moveTowards
// LINESTRING COORDINATES KML FILE
var routePoints  = <?php echo $coordstring;?>
var polycolor    = <?php echo $race_polyclr;?>
var stroke_w = 5;
var polyline = L.polyline(routePoints, {color: polycolor, weight: stroke_w});
map.addLayer(polyline)


texts: Python MapServer MapScript & Mongo - Use GeoJSON String Vs File

texts:
I have been using Python MapScript bindings for MapServer - and they work great. A new datasource for the application I am working on is MongoDB.  I have seen mentions of a MongoDB Driver for MapServer - but that won't work, as some 'data massaging' is required, that is not possible directly in Mongo.

My workflow is this:

Using pymongo - get data of interest from Mongo
'Massage' data - save as valid geojson into tempfile
Mapfile layers are of type MS_OGR - and connection is set to temp geojson file from step above

Overall, this workflow works good, and is stable.  But, some of the geojson payloads are fairly big - say 75-100 MB.  And in some cases, writing geojson to the tempfiles can take a few seconds.  And then I also have to deal with cleanup...

Rather that write the geojson to tempfiles, I am hoping to just use geojson string in the mapfile layers.  But I can't get this work.

Is this even possible?  Any other suggestions?


texts:
I've been trying to install GDAL for a month now in PyCharm. I use PyCharm, I like it, and I'm not switching to/adding Anaconda.

I keep getting the following error (screenshot below) on any GDAL or osgeo-related install. My pip is there and up to date, v. 18.1. Screenshot below of confirming this in command line.


texts:
I have a relatively old laptop (i5 8Gb) but it is fine for most of my GIS needs. Infrequently I would like to run processes which requires larger amounts of RAM (16G+) but I cannot justify forking out lots of money for a high power computer at the moment.

Is it feasible to run software such as GRASS 7.0 or Whitebox-GAT on a rented remote computer/VPS?


texts: Counting how many times .replace() is called during script in python?

texts:
I'm using this code to work on files related to GIS:

import os

path = r"P:/test test"
under_score = [' ', '!', '#', '$', '%', '&', '_-_']
dot = ['.',',']

for path, dirs, files in os.walk(path):
for f in files:
    copy_f = f
    for char in copy_f:
        if (char in under_score): copy_f = copy_f.replace(char, '_')
        if (char in dot): copy_f = copy_f.replace(char, "")
    os.rename(f,copy_f)

Is there any way for me to add a .count() to find out how many characters were replaced when I run this script on a particular path?


texts:
A thing bugging me at the moment is the lack of information on the possibility of finishing a point cloud based 3D model in a desktop software(e.g. Solidworks, Geomagic Design), exporting it in a measurable format (i.e. no meshes or point clouds, but solid bodies, so .sldprt is basically the only one I can export from Solidworks with the textures automatically included), and combining it with map information in WebGIS. For the establishment of my WebGIS, I plan on going the open source way, so a PostgreSQL database and GeoServer are in the picture.

Is there a possibility of including a finished solid model into a WebGIS using open source software?

I haven't found a way which enables me to combine my model with a location, e.g. using OpenStreetMap or a specific WMS service. I know Google Earth has this possibility by uploading a model to the 3D Warehouse, but I can't interact with the model the way I want to there, e.g. measure it. Of course I can do all this in 2D, divide the data into layers and overlap it with map information, but I am really interested in doing it in 3D. To get a more specific idea, I have a gas station model made based on a point cloud, with a lot of pipes, boxes, tiles, stairs,  etc., which I would like to place to its real spatial position.


texts: Using PHP to create a GIS Web application

texts:
I've developed many GIS Applications using ESRI ArcGIS 9.3 and 10.0 for ASP.NET.

I'm wondering if I can do a GIS application using my shapefiles and pure PHP scripting. I've Googled for that purpose, found many MapServer solutions, but they were a little bit complicated.

My question is: can I have a little Web application to mange the shapefiles, create tiles for them, deal with aerial images, and then display map layers using OpenLayers for example?

I need a lite solution (No PHP extensions installing required - NO special set up...etc.).


texts: Only union/dissolve intersecting or adjacent features to speed up query

texts:
I have a table with about 300k polygons, some are overlapping or adjacent to each other and some are islands of their own. I need to dissolve the overlapping/adjacent ones and also keep the non-overlapping ones. I only need the geometries, no attributes.  I can do this with:

create table polygons123_dissolved as
(select (st_dump(st_union(polygons123.wkb_geometry))).geom as geom
from polygons123)

Example:

But this is very slow, more than an hour. How can i speed it up?

Explain output:

Result  (cost=37684.32..37949.59 rows=1000 width=32)
  ->  ProjectSet  (cost=37684.32..37689.59 rows=1000 width=32)
        ->  Aggregate  (cost=37684.32..37684.33 rows=1 width=32)
              ->  Seq Scan on polygons123  (cost=0.00..36538.54 rows=229155 width=849)
                    Filter: ((group_id = 1) AND (type_id = ANY ('{1,2,3}'::integer[])))


texts: How to configure GeoServer user/roles using programmatic way?

texts:
I am working on GeoServer 2.15 on CentOS and configuring users/role manually.
We have 300 layers,10 workspaces & 25 users, created 3 roles (internal, external & classified). I want to assign each role to every workspace and its associated layers.

Manually it's possible, but I am looking at a programmatic way to do this thing.

I am aware its possible using REST interface and curl commands (ultimately its updating that xml file only), so is this a proper way ? Or do I have to write Python code using gsconfig for this thing ?


texts:
I am trying to get a Geoserver wms request using reflect option as static layer to embed into a webpage through a url like this:

http://localhost:8080/geoserver/container/wms/reflect?REQUEST=GetMap&FORMAT=image/png&VERSION=1.1.0&WIDTH=700&HEIGHT=550&SERVICE=WMS&TRANSPARENT=true&LAYERS=layer1&BBOX=Long1,Lat1,Long2,Lat2&SRS=EPSG:4326

,my question is about the zoom level of this layer.
Could be possible to apply a zoom level into this url request?.I would like to apply a zoom level 17.


texts: Arcpy Map Algebra: Maximum Number of Nestings

texts:
When using Map Algebra you can nest multiple expressions as can be seen in the Esri example below.  Is there a limit on how many long of a statement you can have OR a limit on how many Con statements you can nest?  If there is no limit, is there a best practice?

OutRas = Con(InRas1 > 23, 5, Con(InRas1 > 20, 12, Con((InRas1 > 2) & (InRas1 < 17), Sin(InRas1), 100)))


texts: Default resample method used by GEE when exporting images?

texts:
I want to know what the default resample method will be used when I export image like this :

Export.image.toDrive({
  image: image,
  description: 'imagetest',
  **scale: 1000,**
  region:aoi,
  folder: test',
  crs: 'EPSG:4326',
  maxPixels: 9999999999999
});

Nearest？ mean？ cubic?...


texts: Mapserver - Standard Pointsymbolizer edited by SLD

texts:
I have pointfeatures in a OpenLayers map. Every feature has a different attribute or several have even more than one attribute. Thus I would like to use not the standard point symbols like that created by Pointsymbolizer.
Rather I would like to create new points. Is it somehow possible to edit the standard points in the way like filling one point with two colors (one half blue and the other red) by SLD? Or even filling a point with three or four colors?


texts:
Is there a way in QGIS to use an expression to define a polygon feature's color based on whether it's adjacent to other features? I would like to achieve something like the following:

-If a feature has no adjacent features, it gets the default color.
-If two features are adjacent, one gets the default color and the other gets a secondary color.
-If another feature is adjacent to the both of the first two, it gets a tertiary color, and so on.

This is similar to the four-color theorem, except that I want to define the order in which colors are applied, so that the default color is always the first choice and the other colors are given second, third, fourth priority, etc.

The easiest way I can think of is to add a "Color ID" field and manually assign each feature a value of 0, 1, 2 or 3. But can I do this without adding an attribute, just using the existing geometry and searching for shared boundaries?


texts: Finding shared segments between two shapes (line) using QGIS

texts:
I have two line shapefiles. I want to find between them, where are located all shared segments and extract the new selection shapefile.

I am using QGIS. I know that in ArcGIS, there are a tool, named "share line segment with".

Is there some tool like this in QGIS?


texts:
I have multiple lines (lines) intersecting points (centroids) from which I want to extract landuse values and then write the extracted values to the lines as 'land1' and 'land2' values. I need to iterate through the line feature class in order to transfer the appropriate landuse values. My line
layer is made up of 10's thousands of lines so the problem I have is that the select by location operation runs out of memory. I'm using ArcMap 10.0. I have tried processing the dataset in chunks but the problem persists. Why does my 'chunk' approach not work and how can I manage the memory use to avoid this sort of error?

Here is my code:

import arcpy

myfield = "LANDUSE"
obsList = []
obs = arcpy.SearchCursor("lines")
for row in obs:
    myid = row.OBJECTID
    obsList.append(myid)
del row
del obs

chunk = []
mycount = 0
for x in range(0,len(obsList),1000):
    mycount = mycount+1
    chunk = obsList[x:x+1000]
    chunkmax = chunk[len(chunk)-1]
    chunkmin = chunk[0]
    WhereClause = """"OBJECTID">={0} and "OBJECTID"<={1}""".format(chunkmin,chunkmax)
    sle = arcpy.SelectLayerByAttribute_management("lines", "NEW_SELECTION",WhereClause)
    cur = arcpy.UpdateCursor(sle)
    for row in cur:
        counter = 0
        geometry = row.getValue("SHAPE")
        fcs = arcpy.SelectLayerByLocation_management("centroids","INTERSECT",geometry)
        for feat in arcpy.SearchCursor(fcs):
            outval = str(feat.getValue(myfield))
            if counter == 0:
                row.LAND1 = outval
            upcur.updateRow(row)
            if counter == 1:
                row.LAND2 = outval
            upcur.updateRow(row)
            counter = counter +1
        del row
    del cur
    del sle


texts:
I am having difficulty in efficiently dealing with OS Vectormap data that comes in GML format.
In the past I used the Plugin for QGIS (OS GML Tools) to batch convert folders of compressed GML (.GZ) to shapefiles which are then merged and added to the QGIS as a layer. As I now have a lot more data (~150 GML files) and using Shapefiles directly loaded into QGIS is a very slow and inefficient way of displaying/editing the data I have decided to get the data loaded into a PostGIS database.

The first issue is that the plugin doesn't seem to work. When I choose an input folder with only .GZ files (For Example - Vectormap/SP) and then pick the output as an empty folder I created (Vectormap/SP Shapefiles) the process runs and completes but fills the folder with shapefiles only a 1kb in size. When loaded into QGIS these shapefiles appear to do nothing. I have uninstalled the plugin and re-installed to no avail.

If anyone has a better or alternative suggestion for getting GML files into a PostGIS database then I would appreciate it. Surely there musty be an easier way for people using Open Source GIS software to load in Ordnance Survey Vector data?

Thanks for your time.


texts: QGIS 3x viewshed plugin returns wrong results for visibility map

texts:
I am wondering, where I could make the mistake seeing the wrong result, which is coming from the Viewshed plugin in QGIS.
There are a few options, I selected the "Viewshed", which looks like this:

The map generated doesn't match with the same visibility map generated by Heywhatsthat.com service:

whereas the altitude of the point is 172 m.a.s.l. My data attribute table looks as follows:

and accordingly the .csv file, from where my data comes from:

I put height as a m.a.s.l. value, observer height as a 2 meters above ground, radius - 5km.

To be honest I don't understand from where this result comes from.

I found some answer here:
Interpreting viewshed results attribute table

but it doesn't solve my problem, as it refers to the intervisibility analysis.

Can anyone help?

EDITED:

After the brief DEM conversion done here:
Raster -> Conversion -> Translate (Convert format)

I am still getting the wrong result, as per below:

I found a similar problem here:

How to reproject a raster in QGIS?

and I tried this, however the major problem is the target CRS selection:

This is admittedly optional, although without any choice I am getting exactly the same result as above.

After selection the EPSG 27700 OSGB British National Grid 1936 or alternatively another one from Great Britain, I am getting nothing after the Viewshed processing.


texts:
I have identified a need to change VRTs raster order through its order of appearance in it's respective XML file. Unfortunately this seems to not be possible to do dynamically using QGIS commands or plugins.

So I'm thinking of developing a tool that detects a change to the map canvas then gets the center coordinate, finds the rasters in the VRT that are in view then reorders them in the XML file according to how close they are to the canvas center.

This seems very simple to do but I have no idea how to do eventing in the QGIS interface. Ideally this needs to be a background job that continually runs.

Are there any tutorials out there introducing one to eventing in QGIS Python?


texts: Avoid NoData after Focal Statistics Majority

texts:
I am currently working with some classified land cover rasters (containing 8 land cover classes) in ArcGIS 10.0. I would like to make them look smoother and have therefore used the Focal Statistics Majority tool. This results in a raster which, at a first glance, looks great but when I look closer it is filled with thousands of NoData pixels at the boarders between the different land cover classes. I assume this happens because a cell will receive a NoData-value when there is more than one majority value within the neighbourhood. But I wonder if there is any way to avoid this by, for example, making it keep the original value instead of replacing it with NoData when it encounters the problem of having no majority value. Or is there any way to replace the NoData value with the original value after Focal Statistics? Or are there any other recognized, commonly used methods to replace missing data?

Would be very grateful for help!


texts: How to delete and reconnect lines in a polygon in QGIS?

texts:
I made a mistake while I was experimenting clipping and merging polygons and now I have lines that I do not want in my polygon. How can I delete lines which make up the polygon and then connect the gaps?

I would like to get rid of the three lines within this polygon and reconnect the lines at each corner of the triangle.


texts: Apply a mask (shapefile) over a satellite image in Google Earth Engine

texts:
I have a Sentinel 2 satellite image and am rating it. I already classified all the native vegetation separated into a shapefile. Now I need to classify the anthropic uses, I would like the GEE to classify this image but ignore the vegetation areas indicated in the shapefile, how could I do this. I tried this solution Using FeatureCollection to mask an image in Google Earth Engine?, bbut I couldn't replicate,because GEE says "ifc is not a function". How to make it?
Ny code:
var img1 = ee.ImageCollection("COPERNICUS/S2_SR")
.filterMetadata('MGRS_TILE', 'equals', '21KYT')
.filterMetadata('CLOUDY_PIXEL_PERCENTAGE','less_than', 1)
//.filterMetadata('CLOUD_COVER', 'less_than', 0.01)
.filterDate('2020-01-01','2020-12-26')
.select(['B4', 'B3', 'B2'])
.median();

var img2 = img1.clip(table)

//add a região de interesse
var regiao = img1;

//add info iniciais para o treinamento
var treinamento = img1.sample({
 region: table,
 scale: 10,
 numPixels: 450,
});

//var agrupamento = ee.Clusterer.wekaLVQ(9, 0.1).train(treinamento);//6 6,1
//var agrupamento2 = ee.Clusterer.wekaLVQ(9, 0.7).train(treinamento);//6 6,1
var agrupamento3 = ee.Clusterer.wekaLVQ(9, 1).train(treinamento);//6 6,1

//var resultado = img1.cluster(agrupamento);
//var resultado2 = img1.cluster(agrupamento2);
var resultado3 = img1.cluster(agrupamento3);

//Map.addLayer(resultado, {min:0, max: 6, palette: ['#8668d7', '#03a2b8', '#607f72', 'red', '#5d8f5b', '#e624b6', ]}, 'Resultado_Unsup');
//Map.addLayer(resultado2, {min:0, max: 6, palette: ['#8668d7', '#03a2b8', '#607f72', 'red', '#5d8f5b', '#e624b6', ]}, 'Resultado_Unsup2');
Map.addLayer(resultado3, {min:0, max: 6, palette: ['#8668d7', '#03a2b8', '#607f72', 'red', '#5d8f5b', '#e624b6', ]}, 'Resultado_Unsup3');

//Map.addLayer(image, {bands: ['B4', 'B3', 'B2'], min:178, max:1110},  'S2');
// Export the image, specifying scale and region.
Export.image.toDrive({
 image: resultado3.float(),
  description: 'unsup_Sentinel_9',
  scale: 10,
  region: table,
  folder: 'GEE_Ceroula_3',
  //crs: 'EPSG:32721'
  });

var styling = {color: 'yellow', fillColor: '00000000'};
Map.addLayer(table.style(styling));

Map.centerObject(table, 10);


texts: Comments Layer not working in Parks Finder?

texts:
I have installed and begun testing the ParksFinder from ESRI. I have managed to get it running, somewhat, but I keep hitting these weird walls.

My ArcGIS Server is a 10.1 Server.
My SDE DB is a 10.1 Geodatabase on SQL Server.
My WebAdapter is a 10.1 Adapter.
The Parks App I am using is 10.2.

in my config.js on line 132 [URL of CommentsLayer] I am pointing to:
http://maps.cityoftulsa.org/gis/rest/services/LGDM/Parks/MapServer/1

If I change this to a FeatureServer/1 (like the documentation says to) I get an error message whenever I click on a Park:

"Required configuration key values are either null or not exactly matching with layer attributes, This message may appear multiple times."

After this message I get a spinner of doom... the page locks up.

But, if I keep that url pointing to:
http://maps.cityoftulsa.org/gis/rest/services/LGDM/Parks/MapServer/1

Then I can click on a park and in the popup window I can click on the Comments icon, type in a comment, give it a Rating, and click Submit whereupon I get the following error:
"Unable to add comments. Comments table is either absent or does not have write access."

Anyone got any ideas?


texts:
I developed a plugin in QGIS. One of the processes consists in remove some layers from the QGisProject, use them with other third party code in C# (maybe layers are modified by C# code) and finally open again these layers in QGIS. With QGIS 2.18 it works perfectly, but with version 3, when I try to do that, the C# code can't write on removed layers because DBF files are in use by another application (QGIS, I suppose).

I use this code to open layers:

try: #QGis 3.x
       QgsProject.instance().addMapLayer(vlayer, group is None)
except: #QGis 2.x
       QgsMapLayerRegistry.instance().addMapLayer(vlayer, group is None)

Then, to remove layers I use this one:

try: #QGis 3.x
      QgsProject.instance().removeMapLayer(layer.id())
except: #QGis 2.x
      QgsMapLayerRegistry.instance().removeMapLayers([layer.id()])

I think the problem is not in these lines, because some times the process works fine. Other times it works with some layers but not with another layers... Depending of the time that I wait to run the process after open the layers.

If I manually delete the layers from the QGIS user interface and run my process, it works fine.

I have the option of saving the project in a temporal qgis project file. Clear the QGISProject and finally (after C# code) open this temporal file... but it is a dirty solution when I only need to remove one layer...

I don't know what is happening. Some idea? It's possible to free memory of the QGIS project of these removed layers?


texts:
I'm trying to create a buffer as part of a code but I keep getting the following error:

  ExecuteError: Failed to execute. Parameters are not valid. ERROR
  000735: Distance [value or field]: Value is required Failed to execute
  (Buffer).

I set the buffer distance in the script toolbox in ArcGIS but it doesn't seem to recognise it. Here is the code I'm using:

in_features = "Thames"
out_feature_class = "Flood_R_A"
# Run the Buffer Tool
arcpy.Buffer_analysis (in_features, out_feature_class, arcpy.GetParameterAsText(5), "FULL", "ROUND", "ALL", "")

Does anybody know what I'm doing wrong?


texts:
I am generating the SVG file and then use it in jasper report print. When I try to render the map with no sld everything works great, but when I use sld there is this error:

2017-11-30T13:28:36.634+0100  SEVERE  null
java.lang.NullPointerException
    at org.apache.batik.svggen.ImageHandlerBase64Encoder.encodeImage(Unknown Source)
    at org.apache.batik.svggen.ImageHandlerBase64Encoder.handleHREF(Unknown Source)
    at org.apache.batik.svggen.DefaultImageHandler.handleImage(Unknown Source)
    at org.apache.batik.svggen.SimpleImageHandler.handleImage(Unknown Source)
    at org.apache.batik.svggen.SVGGraphics2D.drawRenderedImage(Unknown Source)
    at org.geotools.renderer.lite.StyledShapePainter.renderImage(StyledShapePainter.java:664)
    at org.geotools.renderer.lite.StyledShapePainter.paint(StyledShapePainter.java:226)
    at org.geotools.renderer.lite.StreamingRenderer$PaintShapeRequest.execute(StreamingRenderer.java:3245)
    at org.geotools.renderer.lite.StreamingRenderer$PainterThread.run(StreamingRenderer.java:3572)
    at java.util.concurrent.Executors$RunnableAdapter.call(Executors.java:511)
    at java.util.concurrent.FutureTask.run(FutureTask.java:266)
    at java.util.concurrent.ThreadPoolExecutor.runWorker(ThreadPoolExecutor.java:1142)
    at java.util.concurrent.ThreadPoolExecutor$Worker.run(ThreadPoolExecutor.java:617)
    at java.lang.Thread.run(Thread.java:748)

To generate an SVG I am using https://github.com/geotools/geotools/blob/master/docs/src/main/java/org/geotools/render/GenerateSVG.java (same as here) with added transformation for correct orientation of the image. I use an sld with wkt strings and SVG files. For embedding generated SVG to jasper reports, I had to change the version of batik to 1.8 and to update SVGGraphicFactory to use that version. The problem occurs when generating svg, so I can't convert it to png for printing.

Version of geotools 17.1

EDIT

sample data

https://www.sendspace.com/file/xku43p


texts: How to select behind feature which is overlapped with the same layer on QGIS 1.8.0

texts:
My issue: I received data which was created without "Avoid intersection of new polygon". There a few overlapped feature with the same layer. I could fine overlapped feature by the spatial query. But then, I could not select behind small feature.

At present, I select target feature from attribute table and delete that feature.
I'm hoping I could be able to select feature without attribute table.


texts: Geodjango: Using Lat / Lon to Buffer from Point

texts:
I want to have a radius based distance search. To do this, I want to create a buffer around a point object so that I can filter objects that are inside this buffer, a method that I learned based on a google groups post.

I want to use latitude and longitude in my model definition so that I can easily put the information into my google map later:

class Thing:
    latitude = models.FloatField()
    longitude = models.FloatField()

The only problem is trying to filter results based on lat / lon coordinates:

>>> lat = 37.7762179974
>>> lon = -122.411562492
>>> from django.contrib.gis.geos import Point
>>> pnt = Point(lat, lon)
>>> buf = pnt.buffer(0.0001)

How can I do something like this?

qs = Thing.objects.filter((lat, long)__intersects=buf)


texts: What is Geonode ? What all can be achieved using it?

texts:
Ok i have searched the site and did not find enough information concerning regarding what is a Geonode so i really thought i may ask everyone who has used it to give their opinion.

I am trying to set a website with geospatial data and want to decide whether it is a good solution.

Is it a web application framework like Geodjango, Is it a CMS like drupal for Geospatial information ? Can it be used standalone to deliver both Geospatial and some non geospatial content ?

And honestly , why isn't is as popular as it ought to be ( is it because it is new ?)

The official documentation at http://geonode.org/ provides considerable documentation and examples but i want personal opinions of experts who understand it well.


texts: Automate moving mouse to series of co-ordinates and click

texts:
I have 4 .tif files that I have loaded in QGIS and I am looking to extract data at very specific locations using mutant. I have the points I want loaded in as a csv file and there is a large number of them. Whilst I could manually go to each point and use mutant to locate the data I'm looking for the process can be quite tedious and can be open to errors. Is there a way to automate this process?

Steps I want to automate:

1. Send cursor to specific location
2. Use a click event or manipulate mutant to retrieve data at that point
3. Export that data to a csv or dataframe like structure
4. Iterate over to next location

I have read up on this link which only seems to provide me half the answer and not how to send a click event or manipulate mutant. Is there a way to take control of the mouse pointer using python in QGIS?

I am familiar with python and that is what I intend to use to automate this process


texts: Summing data from a smaller boundary shapefile to a larger boundary shapefile if the boundaries do not align in ArcGIS Pro?

texts:
I have a census tract boundary shapefile for the province of BC, with around 620 polygons, varying in size.

I have another boundary shapefile, BC service areas, this one has much larger polygons and only 28 of them across BC, which do not align with any of the boundaries of the census tract polygons.

I have a population data per census tract ID .DBase file (relating to the census tract polygons), and I have joined this to the census tract boundary file to get the population data per census tract on the map.

However, the end result is that I want to get the census data summarized PER SERVICE AREA. i.e., sum the census tracts within each of the 28 service areas to get the total population of all 28 areas.

I have tried the Summarize Within Tool, this keeps failing, I presume due to boundary overlap errors...

I have thought of converting the census tracts to points or a raster then summarizing that, but I don't think that will work due to duplicating the population data.

Is there a tool that aligns the Service Area Boundary layer to the census tract boundaries? Or an overall process to get the census data summarized per service area?


texts: Way processing is too slow in osm2pgsql?

texts:
I'm importing the planet data file. After see the Node processing taking about a hour at 2800.0 k/s to finish I think all job will be fast but now Way is still working after 28 hours at 0.5 k/s and done just 31.583k of 402.887k ! (OSM stats).

The top command shows it taking 15% of RAM and 5% of CPU.

My import command was

$osm2pgsql --number-processes 6 --flat-nodes ./osm_flat_nodes.db --disable-parallel-indexing --latlong --verbose --create --hstore --slim --cache 4000 --database osm --username postgres --host 127.0.0.1 --style ./default.style ./planet-latest.osm.pbf

My machine:
CentOS / 32GB RAM / 600GB HDD / 8 Cores

Is this normal?


texts:
I have some data in a raster stack (CCSM4 climate projections from 2010 - 2300) that I've done some processing on, but now I would like to change the extent. I'm not 100% sure how to do this in R. So far I've tried Extent, Origin, projectExtent, projectRaster, but haven't had any luck. Below are two images showing the issue I'm having and the results of what I've tried, and a reproducible example

Here is my original dataset, the extent is

> extent(pr_ccsm4_2.5deg)
class       : Extent
xmin        : -0.625
xmax        : 359.375
ymin        : -89.5288
ymax        : 90.4712

What I would like to do is get the data to have an extent of

> extent(ras25)
class       : Extent
xmin        : -180
xmax        : 180
ymin        : -90
ymax        : 90

I thought that I could simply change the extent using extent(pr_ccsm4_2.5deg) <- extent(ras25), however this results in my data looking like this

So now, while the extent is correct, the data is not projected properly.

Here is a reproducible example

library(rworldmap)
library(cleangeo)
library(raster)
library(rasterVis)

# continental shapefile
cont <- getMap()
cont <- clgeo_Clean(cont)
cont <- sapply(levels(cont$continent),
               FUN = function(i){
                 poly <- gUnionCascaded(subset(cont, continent == i))
                 poly <- spChFIDs(poly, i)
                 SpatialPolygonsDataFrame(poly, data.frame(continent = i, row.names = i))
               }, USE.NAMES = TRUE)
cont <- Reduce(spRbind, cont)

ras <- raster(xmn = 0, xmx = 360, ymn = -90, ymx = 90, res = 2.5)
values(ras) <- runif(n = ncell(ras), 0, 1)
ras
> ras
class       : RasterLayer
dimensions  : 72, 144, 10368  (nrow, ncol, ncell)
resolution  : 2.5, 2.5  (x, y)
extent      : 0, 360, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : layer
values      : 2.22486e-05, 0.9998961  (min, max)
# Wrong extent
rasterVis::levelplot(ras, margin = NA, par.settings = RdBuTheme) + layer(sp.polygons(cont))

ras2 <- raster(xmn = -180, xmx = 180, ymn = -90, ymx = 90, res = 2.5)
values(ras2) <- runif(n = ncell(ras2), 0, 1)
ras2
> ras2
class       : RasterLayer
dimensions  : 72, 144, 10368  (nrow, ncol, ncell)
resolution  : 2.5, 2.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : layer
values      : 0.0001686567, 0.9998712  (min, max)

# Correct extent
rasterVis::levelplot(ras2, margin = NA, par.settings = RdBuTheme) + layer(sp.polygons(cont))

rasPrj <- projectRaster(ras, ras2, over = FALSE)
rasPrj
# Correct extent, but missing data
rasterVis::levelplot(rasPrj, margin = NA, par.settings = RdBuTheme) + layer(sp.polygons(cont))


texts: Way of viewing attribute table of layer at end of model without loading layer in QGIS

texts:
I have created a model in QGIS 3 graphical modeler, with the final output being a shapefile. I was wondering if it is possible to have the final result as a pop up of the data table just to view the result of the model. Ideally I don't even want the shapefile, even though you can load it as a temporary layer.

The final attribute table looks like this:


texts: Merging multiple raster files incorrect value in ArcGIS

texts:
I have multiple rasters in .tif (geotiff) format and I´m trying to merge them into one using the Mosaic to New raster tool in ArcGIS(10.7.1). These rasters have a partial overlap between them and are 16bit unsigned int values.
When i merge them using the Mosaic to New raster tool with proper folder locations, the merged file name (something.tif), pixel type (16 bit unsigned int), number of bands (4 - R, G, B and NIR), Mosaic operator (last) and mosaic colourmap (first) selected i get a new raster whose values are sometimes different than the underlying rasters in a few pixels.
For example, these are the values for a colocated pixel:
Merged raster: 919,1182,619
Raster 1: 883,1124,631
Raster 2:912,1163,668
The merged raster values are very close to the underlying raster in some pixels and same values in some pixels but shouldn´t it be the same as the last raster encountered throughout the image? Why is it different?
I´ve attached images of these pixel values along with some info on the merge function i ran.


texts: Is there any accuracy difference between standard DTM editing tools like DT Master and automatic removal of artefacts using ERDAS?

texts: I would like to know whether there is any accuracy difference while editing DTM points using standard photogrammetry software like Inpho DT Master and automatics removal of artefacts in ERDAS imagine.


texts: QGIS processing load memory layer with specific name

texts:
I have a simple QGIS plugin and I use processing algorithms from my models.

Now I want to load my export from processing in my QGIS project that I do easy with runandload but I can't to define name that show me an default name how to load memory layer with my name?

outputs_QGISFIELDCALCULATOR_4=processing.runandload('qgis:fieldcalculator', outputs_QGISFIELDCALCULATOR_3['OUTPUT_LAYER'],'newnew',2,50.0,0.0,True,'left(\"TOP_KOIN\", strpos(\"TOP_KOIN\" ,\' \')-1)',"memeory:layer")

And if I need that memory layer to continue my script for example to do something else how to take ?with memory name? or like this outputs_QGISFIELDCALCULATOR_4?

update

if result:
    selectedLayerIndex = self.dlg.comboBox.currentIndex()
    selectedLayer = layers[selectedLayerIndex]
    print selectedLayer
    #calc='C:/Users/saran/Desktop/qq/calc.shp'
    outputs_QGISFIELDCALCULATOR_1=processing.runalg('qgis:fieldcalculator', selectedLayer,'trid',1,10.0,0.0,True,'$id',None)
    load_it1 = QgsVectorLayer(outputs_QGISFIELDCALCULATOR_1['OUTPUT_LAYER'], 'tid', 'ogr')
    outputs_QGISFIELDCALCULATOR_2=processing.runalg('qgis:fieldcalculator', outputs_QGISFIELDCALCULATOR_1['OUTPUT_LAYER'],'area',0,10.0,3.0,True,'$area',None)
    load_it2 = QgsVectorLayer(outputs_QGISFIELDCALCULATOR_2['OUTPUT_LAYER'], 'area', 'ogr')
    outputs_QGISFIELDCALCULATOR_3=processing.runalg('qgis:fieldcalculator', outputs_QGISFIELDCALCULATOR_2['OUTPUT_LAYER'],'perimetre',0,10.0,3.0,True,'$perimeter',None)
    load_it3 = QgsVectorLayer(outputs_QGISFIELDCALCULATOR_3['OUTPUT_LAYER'], 'perimetre', 'ogr')
    outputs_QGISFIELDCALCULATOR_4=processing.runandload('qgis:fieldcalculator', outputs_QGISFIELDCALCULATOR_3['OUTPUT_LAYER'],'newnew',2,50.0,0.0,True,'left(\"TOP_KOIN\", strpos(\"TOP_KOIN\" ,\' \')-1)',None)

    load_it4 = QgsVectorLayer(outputs_QGISFIELDCALCULATOR_4['OUTPUT_LAYER'], 'new', 'ogr')
    QgsMapLayerRegistry.instance().addMapLayer(load_it1)
    QgsMapLayerRegistry.instance().addMapLayer(load_it2)
    QgsMapLayerRegistry.instance().addMapLayer(load_it3)
    QgsMapLayerRegistry.instance().addMapLayer(load_it4)
    my_output1 = processing.getObject(outputs_QGISFIELDCALCULATOR_1['OUTPUT_LAYER'])
    my_output2 = processing.getObject(outputs_QGISFIELDCALCULATOR_2['OUTPUT_LAYER'])
    my_output3 = processing.getObject(outputs_QGISFIELDCALCULATOR_3['OUTPUT_LAYER'])
    my_output4 = processing.getObject(outputs_QGISFIELDCALCULATOR_4['OUTPUT_LAYER'])


texts: How to create "undulating terrain" for playground in ArcGIS?

texts: How could I create undulating terrain in ArcGIS? I know I could create it by creating a TIN model but first I need to have "input data" - shapefile with Z (height in attribute). Is there any way I can create something you see at the picture below? The shapefile will be used in ArcScene for 3D of course.


texts: Blank map in result page using MapServer

texts:
I am using MapServer to retrieve the data into my map but I get only a blank page instead of data on the map.

I am using the below query to fetch the data from a PostgreSQL:

DATA "the_geom from aoi using unique gid using SRID=2263"

Is this the correct way to form the query in the Mapfile?


texts: Copy a geometry from one feature to another without attributes?

texts:
With QGIS it's easy to add a new part to an existing feature/geometry but I can't see how to copy/paste this part from another existing feature/geometry.

Is it possible?


texts: Wrong scale units when reading raster grib2 in QGIS 3

texts:
I downloaded grid2 files to read meteorogical variables.
I wanted to read them through QGIS 3.10.2 (windows and I have to stick with this version) but I encountered problems with the scale of the variables. It should be between 0 and around 150, but I have values around 0 and 560,000. There are missing values that are coded as 9999 but those missing values seems to be correctly accounted for when looking into the QGIS metadata of the grib2 files.
I do not have this problem when using QGIS 2.18, and when reading the grib2 file through python (using pygrib package) I still have no problems (values are in the right magnitude).
It seems to be linked to the way I read the file in QGIS 3.