texts: WFS layer not printed in OpenLayers Map?

texts:
This is the WFS service I'm using:

  https://idena.navarra.es/ogc/wfs?service=WFS&version=1.1.0&request=GetFeature&typename=IDENA:DOTACI_Sym_AlojTur&outputFormat=application/json&srs=EPSG:25830

I add it to my OpenLayers map as a Vector layer as following:

proj4.defs("EPSG:25830","+proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs");
register(proj4);
const proyeccion = getProjection('EPSG:25830');

var vectorSource = new VectorSource({
  format: new GeoJSON(),
  url: function(extent) {
    return 'https://idena.navarra.es/ogc/wfs?service=WFS&version=1.1.0&request=GetFeature&typename=IDENA:DOTACI_Sym_AlojTur&outputFormat=application/json&srs=EPSG:25830';
  },
  strategy: bboxStrategy
});

var vector = new VectorLayer({
  source: vectorSource,
  projection: 'EPSG:25830',
  style: new Style({
    stroke: new Stroke({
      color: 'rgba(215, 29, 29, 1.0)',
      width: 4
    })
  })
});

var map = new Map({
  target: 'map',
  layers: [vector,layer], //layer is an WMS layer
  view: new View({
    projection: proyeccion,
    center: fromLonLat([-1.629950,42.63],proyeccion),
    zoom: 9
  })
});

When loading the map, the WFS request is done and the response contains the feature collection in GeoJSON format and with the correct EPSG:

  Although everything seems to work fine, the map does not show the wfs layer

What is the problem here?


texts: Attribute table of shapefile completely blank and empty after a join - ArcMap - how to recover it?

texts:
I'm using ArcMap 10.6. I have a shapefile of geo-political borders from Natural Earth. Both the shapefile and its attribute table are in working order. But after I join a csv to the shapefile, the attribute table of the shapefile completely disappears. I can't figure out where it's gone or how to recover it. There are absolutely no columns or rows (see photo). I know the join worked, because if I click/identify an individual polygon of the shapefile, I can see its attributes (new and old) using the identify tool; I can also see the fields in the properties > fields window, but no matter what I have tried, I cannot see anything in the attribute table of the shapefile after the join. The menus in the attribute table also are not working. However, if I remove the join, the attribute table of the shapefile once again appears as normal.

I tried exporting the data after the join, as another shapefile - to see if the attribute table would be visible there - but the export progress remains at 0% after several minutes...so I think that's not going to work.

How can I make the attribute table visible again after I do the join?


texts: Debugging Python script using Visual Studio in ArcGIS Pro?

texts:
According to these links, one can debug a Python script using Visual Studio in ArcGIS Pro.
Debug Python code—ArcPy Get Started | ArcGIS Desktop
ArcGIS Pro and Python - YouTube
I am using ArcGIS pro version 2.3 and Visual studio 2017.
My steps

Downloaded the 'Python development' workload to visual studio Pro 2017.
Edited my .py file to add some intentional errors for testing.
Opened ArcGIS Pro with my tool and my data to be run.
Opened my .py file in Visual Studio and added some breakpoints.
Attached to process
in Visual Studio. See image below:

Ran my script tool in ArcGIS Pro.

The tool fails in ArcGIS Pro because it has errors but I am not re routed to Visual Studio to see the error line as I expected. I then removed the errors in the code and left only the breakpoints but I still wasn't able to hit the breakpoints. Again, the tool completed without re-routing me back to visual studio. Has anyone been able to successfully debug Python tools using Visual Studio Pro 2017 in ArcGIS Pro?


texts: Keep spatial lines intersecting spatial polygons using R?

texts:
I am trying to keep spatial Lines objects intersecting a spatial polygon in R. The intersect function from raster package is doing actually a crop. But what I need is to keep the lines intersecting a polygon (but not to crop them). In ArcMap there was a possibility for that (something like select by attributes using a polygon).

Some code I created:

library(sp)
library(raster)

# from the sp vignette:
l1 = cbind(c(-35,-32, -31, -18, -15, 0, 22, 24 ), c(38, 39, 41, 44, 45, 48, 55, 66))
rownames(l1) = letters[1:8]
l1a = cbind(l1[,1]+25,l1[,2]+25)
rownames(l1a) = letters[1:8]
l2 <- l1-66
rownames(l2) = letters[1:8]
l3 <- l2 - 28.2
rownames(l3) = letters[1:8]

Sl1 = Line(l1)
Sl2 = Line(l2)
Sl3 = Line(l3)
S1 = Lines(list(Sl1), ID="a")
S2 = Lines(list(Sl2), ID="b")
S3  = Lines(list(Sl3), ID  = 'c')
Sl = SpatialLines(list(S1,S2, S3))

poly <- extent(c(-45, 15, -15, 57))
poly <- as(poly, 'SpatialPolygons')

plot(Sl, col = c("red", "blue", 'green'))
plot(poly, add  = T, lwd = 2) # here I should keep the red and blue lines

intersected <- raster::intersect(Sl, poly)

plot(intersected)
plot(poly, add=T) # I managed only to crop the lines which is not exactly what I wanted.


texts:
Using ArcGIS 10.4, arcpy

I am having trouble understanding field mappings and the esri examples aren't helping at all.  When I use model builder the field mappings make sense, but I need to customize them more.  I can successfully take the ObjectID field,add it to my new feature class with the input data from the old, but that is the only field that comes in. Is there a way to add all the others instead of going one by one? It seems like it should be so simple, I tried the (addtable) but that is where I am lost.

My end game is to keep ObjectID in my new feature class, I want everything else the same (Type, Length, etc).  Right now on my 'FeatureClasstoFeatureClass' all I can keep with the field mappings is ObjectID.  I could join it up and rexport but field mapping seems like it should be the easiest.

http://pro.arcgis.com/en/pro-app/arcpy/classes/fieldmappings.htm

I am looking at the examples here.  They only mention adding two fields, but what about every field.


texts:
The National Weather Service (NWS) makes the results from its Global Forecast System model (GFS) available for download as GRIB2 files.

While accessing current (or recent) forecasts is straight-forward, I am looking for the archived forecasts over at least the last year, better multiple years back for the purpose of statistical analysis of predominant weather patterns in certain areas.

Does anybody know how to access this information?


texts: Can QGIS be used to define territories of equal population?

texts:
I was wondering whether the QGIS plugin Cmaps region generator can be used to define territories of equal population by splitting a file of small postal areas/polygons and then automatically building territories up with an equal size population?

Or are there other open source software that could do this?


texts:
This is a find and replace script I have been working on. I feel I am very close, as the script does everything but actually change the values, leading me to believe there is an operator issue. I have been using == as you can see, but have also tried = when the cursor looks for oldText only to get a syntax error. Any thoughts?

import arcpy

FC = arcpy.GetParameterAsText(0)
targetField = arcpy.GetParameterAsText(1)
oldText = arcpy.GetParameterAsText(2)
replaceText = arcpy.GetParameterAsText(3)

fields = arcpy.ListFields(FC)

queryString = queryString = targetField + ' = ' + oldText
print queryString
arcpy.AddMessage(queryString)

try:

    with arcpy.da.UpdateCursor(FC, (targetField), queryString) as cursor:
        for field in fields:
            if field.name == targetField and field.type == 'String':
                for row in cursor:
                    if row[0] == str(oldText):
                        row[0] = str(replaceText)
                    cursor.updateRow(row)

        elif field.name == targetField and field.type == 'Single':
             for row in cursor:
                if row[0] == float(oldText):
                    row[0] = float(replaceText)
                    cursor.updateRow(row)

        elif field.name == targetField and field.type == 'Integer':
             for row in cursor:
                if row[0] == int(oldText):
                    row[0] = int(replaceText)
                    cursor.updateRow(row)

        arcpy.AddMessage("String where " + oldText + " found and replaced with " + replaceText)

except:
        arcpy.AddMessage \(arcpy.GetMessages())
        arcpy.AddError("Failed to find and replace where " + queryString)


texts: Setting Script Tool Parameters Using Python and ArcGIS 9.3.1?

texts:
Can't quite figure this out and i'm sure it's a simple fix.  I have a script that utilizes the AGSSOM tool.  It first stops all published services, allows for databases to be tidied up by Compressing/Uncompresing and Compacting, then resarts the services.  The script is a tool in a tool box.  I'm currently trying to set a parameter where the user picks which geodatabases to clean up (our services use several geodatabases and they are constantly being manipulated).  Within the Parameters tab of the scripts properties, I've added a parameter and assinged it to Data Element as the Data Type (anyone have a better suggestion?).  Additionally, I've set the MultiValue to Yes.  Here is an example of some of the code.

inGDB = gp.GetParameterAsText(0)
gp.CompressFileGeodatabaseData_management(inGDB)

When I load multiple GDBs It crashes whith an error:

 <class 'arcgisscripting.ExecuteError'>: Failed to execute. Parameters are not valid.
Input file geodatabase data: Dataset 'Z:\Data\Reference.gdb;Z:\Data\As-Built.gdb' does not exist
Failed to execute (CompressFileGeodatabaseData).

Failed to execute (ArcgisServerUpdater).

If I load a single GDB It works.


texts: What data format to do temporary editing with QGIS before committing it to a RDBMS?

texts:
I am trying to create a temporary spatial data from a PostGIS database that a user can do the editing on. Currently, I have a service that returns GeoJSON from the PostGIS database. But I don't think it's possible to edit GeoJSON using QGIS. What is the reason?

Secondly, what would be the best format to do the editing then to convert it back to GeoJSON. I have shapefile and Spatialite in mind. The downside of shapefile is that I can only use one table per one shapefile. If for example I have a line and a polygon table, I would need two shapefiles for each. I am not sure what advantage Spatialite offers over shapefile? Is it possible to edit Spatialite data on QGIS? How comparable is Spatialite to Esri Personal Geodatabase?


texts:
Here is what I do in django:

class MyModel(models.Model):
    pass

class MyModelAdmin(admin.ModelAdmin):
    pass

admin.site.register(MyModel, MyModelAdmin)

But with geodjango I have this:

  admin.site.register(MyModel, admin.OSMGeoAdmin)

So I to create my MyModelAdmin class with geodjango?
I want to add a search_fields


texts: Problem when trying to cut linestring at multiple points after conversion from Multi-linestring

texts:
This is my script:

def cut_line_at_points(line, points):
    '''
    line is gdf format
    points are series of points in gdf format
    '''
    # union all geometries - convert the gdf to shpely (multi)line and (multi)point
    line = line.geometry.unary_union
    points = points.geometry.unary_union

    # check if line in multi-part and merge all parts to one linestring
    if line.geom_type != 'LineString':
        # Put the sub-line coordinates into a list of sublists
        outcoords = [list(i.coords) for i in line]
        # Flatten the list of sublists and use it to make a new line
        line = LineString([i for sublist in outcoords for i in sublist])

    # First coords of line
    coords = list(line.coords)

    # Keep list coords where to cut (cut = 1)
    cuts = [0] * len(coords)
    cuts[0] = 1
    cuts[-1] = 1

    # Add the coords from the points
    coords += [list(p.coords)[0] for p in points]
    cuts += [1] * len(points)

    # Calculate the distance along the line for each point
    dists = [line.project(Point(p)) for p in coords]

    # create new df and sort again if geom_type was LineString

    # sort the coords/cuts based on the distances
    coords = [p for (d, p) in sorted(zip(dists, coords))]
    cuts = [p for (d, p) in sorted(zip(dists, cuts))]

    # generate the Lines
    lines = []
    for i in range(len(coords)-1):
        if cuts[i] == 1:
            # find next element in cuts == 1 starting from index i + 1
            try:
                j = cuts.index(1, i + 1)
                lines.append(LineString(coords[i:j+1]))
            except:
                cuts[-1] = 1
                j = cuts.index(1, i + 1)
                lines.append(LineString(coords[i:j+1]))

    # transform Geometry Collection to GeoDataFrame
    segments = [feature for feature in lines]
    # set the geometry and other columns
    gdf_segments = gpd.GeoDataFrame(list(range(len(segments))), geometry=segments)
    gdf_segments.columns = ['index', 'geometry']

    return gdf_segments

everything is working really well when it's a shapely.LineString

but when I'm trying to cut the Multi-linestring after the flattening it all goes south:

I found out that when Calculating the distance along the line it scrambles the points distances, in this line:

dists = [line.project(Point(p)) for p in coords]

but I don't have any idea how to fix it. I use Python 3.7 and Shapely.


texts: QGIS 3.4 Modeler gives error on non-required parameter?

texts:
I created a model which is asking to input 20 OPTIONAL parameters for vector layers.

But when the user does not enter a value for any of the 20 parameters, QGIS returns an error.

Traceback (most recent call last):
File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\algs\qgis\PointsLayerFromTable.py", line 91, in processAlgorithm
raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT))
_core.QgsProcessingException: Cannot load source layer for INPUT: invalid Value

Note: I translated the last line.

Is there a way, to not enter a value for some of the parameters and still have the ones processed that do have values?


texts: How to represent an imaginary flat world in QGIS without using a round-world CRS

texts:
Like one or two other users, I'm interested in using QGIS to map an imaginary flat world. A previous answer suggested simply refraining from defining a CRS for such a world, but I cannot seem to make that work well with QGIS 1.7 on Mac OS X. Pretty much every time I interact with a layer, the software warns me that no CRS is defined, and if I ever inadvertently open the CRS dialog for any reason, QGIS appears to apply the default CRS to my map. That happens quite often when I'm fumbling around looking for a way to find a setting that will establish the appropriate scale for distances on the map.

The fictional world is a flat plane of unknown (and possibly infinite) extent. The map is centered on a location that is locally taken to be the center of the world. Although the world is of unknown extent, the mapped area is about ten thousand kilometers north-to-south, and the same east-to-west.

I can map the world adequately using tools like Photoshop and Illustrator, but obviously they aren't GIS systems. I'd like to use QGIS for the myriad additional data-handling features it would provide, but I'm pretty well stumped by my inability so far to find or create a CRS --or to get QGIS to accept working without one.

I'm close to concluding that I'll have to write my own quasi-GIS system for this purpose that obeys the (admittedly idiosyncratic) rules of the fictional world.


texts: Creating Voronoi diagrams using JTS and the manhattan or minkowski distance metrics

texts:
In Java I am using the JTS library. Is it possible to use this library to generate Voronoi sites using different types of distance functions?

https://locationtech.github.io/jts/javadoc/org/locationtech/jts/triangulate/VoronoiDiagramBuilder.html

for example
*voronoi sites using Euclidean Distance:

*voronoi sites using Manhattan Distance: ( what i want to do using JTS but dont know how )


texts: How to change a project variable using Python?

texts:
I'd like to change the value of a user-defined project variable (can be manually edited in Project Settings | Variables) using the Python console. I tracked down the setVariable() function in the QgsExpressionContextScope class but haven't succeeded in actually changing the variable in the project settings. My code so far:

iface.mapCanvas().mapSettings().expressionContext().scope(0).setVariable('myvar',1)

I guess I'm getting lost in the different expression contexts ...


texts:
I would like to know the area that each value in a raster has... The raster I have has values from 2-15 and I would like to know the area that each value is occupaying within my polygon... Thanks!

I vectorized and calculated the area using geometry columns but I find a problem when I add all of the calculated segments, they don't match the poligon's total area!

The fist picture is the resulting raster, I clipped my polygone out and then vectorized it (pic of de right)...


texts:
My script: https://code.earthengine.google.com/03072a1aa4dd0c8bcb5f0d6f50840cc7

This function used to return a token. Now token is empty and from what I can remember the mapId used to look different as well.

This is what getMap() gives me:

{
  "mapid": "projects/earthengine-legacy/maps/6549925a95e57c05d203d65bc6820c40-e563a4b747bc6703107afb69bbe17a7e",
  "token": "",
  "image": ...image stuff
}

How should i format my tile url now?


texts: Google engine MODIS and filter by polygon

texts:
I tried to apply "filterBound" function to a MODIS image collection, using a polygon, but it doesn´t work.

This is my script :

var monte = ee.FeatureCollection("ft:1Sm4ObZgrKeCP-VlFa6Jf019-62oDS_AfqvqYw0t6"),

    MOD_W = ee.ImageCollection("MODIS/MOD13Q1");
var MOD = MOD_W.filterDate('2002-07-04', '2017-03-30');

var MOD_monte = MOD.filterBounds(monte);

Map.addLayer (MOD_monte);

When I use the same script for Landsat collection, it works.

How can I do for apply filterbounds to MODIS?


texts: How to auto Correct Polygons that are not drawn correctly?

texts:
I was part of a project recently where handheld GPSs were used to collect vertices of parcels. After plotting the points using 'P2One' plugin in QGIS I found out that some of the polygons were twisted due to the fact that the vertex points were not collected in a clock-wise or counter-clockwise manner. The 'P2One' plugin uses the time stamps for each point on the GPS data collected to determine the sequence with which to plot points collected for a parcel.

Is there a tool or script in QGIS that can help correct the twisted polygons to regular shapes automatically?


texts: Exporting Shapefile to SpatiaLite format using QGIS?

texts:
I am trying to export a shapefile of 3221 mixed POLYGON and MULTIPOLYGON geometry elements and every time, I get this same error:

Does anyone have any idea what's going on?

And no, this is not an exact duplicate of https://stackoverflow.com/questions/15203403 as I am working solely within QGIS.


texts: Refreshing Versioned Feature Class using ArcObjects?

texts:
I have a FeatureClass A.
This Feature Class is Versioned.

User X has a tool open, that reads records from this Feature Class, by clicking a refresh button. The refresh button code is basically a query to the table using this statement

FeatureClassA.Search(null,false);

Another user Z decided to edit the table in a different session. So they open an arcmap session, adds the feature class, starts editing, and creating features.

Once the user finished editing, they save the edit session.

The first user X, clicks the refresh button again, using the same query above, no new Features are brought besides the ones that were already there.

If the user closes the edit session, and reopens it, the new record will show.

How can I make the refresh button get the latest changes from the Feature Class without having to stop and restart the edit session?


texts: How to include long.vrt and lat.vrt into bt.vrt as x-dataset and y-dataset respectively, Using python?

texts:
I am using following code to create bt2.vrt file using python.

import numpy as np
import gdal
from gdalconst import *
import osr

# Function to read the original file's projection:

def GetGeoInfo(FileName):
    SourceDS = gdal.Open(FileName, GA_ReadOnly)
    xsize = SourceDS.RasterXSize
    ysize = SourceDS.RasterYSize
    GeoT = SourceDS.GetGeoTransform()
    Projection = osr.SpatialReference()
    Projection.ImportFromWkt(SourceDS.GetProjectionRef())
    DataType = SourceDS.GetRasterBand(1).DataType
    #DataType = gdal.GetDataTypeName(DataType)
    return xsize, ysize, GeoT, Projection, DataType

# Function to write a new file.

def Create_vrt(Name, Array, driver,xsize, ysize, GeoT, Projection, DataType):
    if DataType == 'UInt16':
        DataType = gdal.GDT_Float32
    NewFileName = 'bt2'+'.vrt'
    # Set nans to the original No Data Value
    # Set up the dataset
    DataSet = driver.Create( NewFileName, xsize, ysize, 1, DataType )
            # the '1' is for band 1.
    DataSet.SetGeoTransform(GeoT)
    DataSet.SetProjection( Projection.ExportToWkt() )
    # Write the array
    DataSet.GetRasterBand(1).WriteArray( Array )
    return NewFileName

# Open the original file
FileName='/home/vigna/TIFFs/bt.tif'

DataSet = gdal.Open(FileName, GA_ReadOnly)
# Get the first (and only) band.
Band = DataSet.GetRasterBand(1)
# Open as an array.
Array = Band.ReadAsArray()

# Now I do some processing on Array, it's pretty complex
NewArray = Array
print(NewArray)

#save the new file, in the meantime I have
xsize, ysize, GeoT, Projection, DataType = GetGeoInfo(FileName)
# Set up the vrt driver
driver = gdal.GetDriverByName('VRT')

# Now turn the array into a vrt.
NewFileName = Create_vrt('bt3', NewArray, driver,
                        xsize, ysize, GeoT, Projection, DataType)
#Create_vrt()
 print("ok")

bt2.vrt file created using above code is : [1]: https://i.stack.imgur.com/JYqVf.png

According to  [1]: Unable to warp HDF5 files
I need to include lat.vrt and long.vrt as y_dataset and X_dataset respectively to perform "gdal_warp" on bt.vrt file.
How to include it using python code?


texts: What might cause ArcGIS Online Public Account to reset and lose all content?

texts:
I have two ArcGIS Online accounts:

ArcGIS for Organizations which I log into as a named user
Public Account which I log into as a Gmail user

Today, I logged into my Public Account and as I did so I saw a message which said something about "resetting" before the login proceeded.

When I visited My Content I found that all of the considerable content that I authored there last week was gone.

If anybody else has experienced this can you tell me the likely/precise wording of the "resetting" message that I saw, so that I can submit this to my local Esri support to look into?

I now think the intermittent issue above is probably unrelated to the reset message which I have since seen again without loss of content:

  Your login session has been reset. Please try again.


texts: Formula for converting SWEREF99 to WGS84?

texts:
I have a few hundred rows of data with sweref99 coordinates that I want to plot using Fusion Tables on Google Drive.

What formula can I use to convert the two columns of Sweref99 x and y coordinates?

All my Google results so far just makes me end up at converters, not any formulas, forcing me to manually convert every row.


texts: How to select multiple Z field from attribute table to perform IDW interpolation using ArcPy

texts:
I have 12 point feature datasets for different month. I want to perform IDW interpolation.

I have created the python code but  its taking only one "Z" field for different month feature datasets from the attribute. Below I have mentioned the code. How to assign here multiple Z field from each month folder (JANUARY, FEBRUARY......DECEMBER).

Point feature class from different month folder (e.g. For January Merged_001, for FEBRUARY Merged_002.......for DECEMBER Merged_012 ) and want to perform IDW interpolation and save it as on same month folder as it Z field name (Max_Temper, Min_Temper, Precipitatat, Wind, Relative_H, Solar) with month  (e.g for January 001_Max_Temper, 001_Min_Temper, 001_Precipitatat, 001_Wind, 001_Relative_H, 001_Solar).

After  Whatever Solar raster for different month  (for January 001_Solar....for DECEMBER 012_Solar) we will get, I want to do some calculation (001_Solar * 30) *0.5 and save it 001_sr ........012_sr in Month wise folder.

Here is my code.

import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.env.parallelProcessingFactor = "100%"
arcpy.env.overwriteOutput = True

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

env.workspace = "D:\SWAT-WEATHER-DATA2\APRIL"

# Set local variables
inPointFeatures = "D:\SWAT-WEATHER-DATA2\APRIL\Merged_004.shp"

zField = "Min_Temper" #Z_Filed Name= Max_Temper, Min_Temper, Precipitatat, Wind, Relative_H, Solar

cellSize = 0.002298707671
power = 2
searchRadius = RadiusVariable(10, 150000)

#Mask region of interest
mask="D:\Gujarta Shape file\GUJARATSTATE.shp"

# Execute IDW
outIDW = Idw(inPointFeatures, zField, cellSize, power, searchRadius)

# Execute Mask
IDWMASk = ExtractByMask(outIDW, mask)

# Save output, except Solar Radiation raster
IDWMASk.save("004_Min_Temp.tif")

# 004_Max_Temp.tif, 004_Min_Temp.tif....Solar

#Only for Solar Radiation raster
#PAR = (IDWMASk * 30) * 0.5

#Save output, Only for Solar Radiation raster
#PAR.save("003_sr.tif")

#print done


texts: Converting double field to integer in database view so ArcMap query layer can use field as unique identifier

texts:
I have a field in a table that is a double data type, but the values are actually integers. The table is from third-party software, so I don't have control of the field data types.

I'm attempting to make a database view of the table, and add it to ArcMap (which I believe makes it into a query layer).  Arcmap requires a unique identifier field be selected when the view is added to the map. There isn't a single unique field in the table, but the combination of two other fields plus the double field is unique. I plan to select all three fields as the unique id (sort of like a composite key). The problem is that it that ArcMap seemingly does not recognize double fields as unique identifiers (and probably for good reason). Text fields and integers seem to be recognized just fine.

As a workaround, I'm trying to convert the double field to an integer field in the database view. I've tried various combinations of the TO_NUMBER, ROUND, FLOOR, TRUNCATE, and CAST Oracle functions, with no success.

How can I get the table into ArcMap as a database view, and use the double field as a unique identifier?

Environment: Oracle 11g, SDE 10.3.1, ArcGIS Desktop 10.3.1

Update Jul-29-2016:

The table is non-spatial. It will just be used to join to another non-spatial table about once per quarter for QC purposes. The table has ~20,000 records. There isn't a ROWID or OBJECTID field in the table, thus forcing me to find a workaround.


texts:
I have downloaded QGIS and loaded the Google Layers plugin.  I have imported a DXF CAD drawing.  Question is how do I place the shapefile in the correct place on the map?

Also, how do I tell Google Layers to return a particular area of the map?


texts:
I was wondering how to implement a directed (or one-sided) Hausdorff distance function in PostGIS?

I looked into the PostGIS 2.3.3 source code, and it seems that ST_HausdorffDistance is implemented as a simple C wrapper function that calls the corresponding GEOS function as follows:

PG_FUNCTION_INFO_V1(hausdorffdistance);
Datum hausdorffdistance(PG_FUNCTION_ARGS)
{
    GSERIALIZED *geom1;
    GSERIALIZED *geom2;
    GEOSGeometry *g1;
    GEOSGeometry *g2;
    double result;
    int retcode;

    POSTGIS_DEBUG(2, "hausdorff_distance called");

    geom1 = PG_GETARG_GSERIALIZED_P(0);
    geom2 = PG_GETARG_GSERIALIZED_P(1);

    if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
        PG_RETURN_NULL();

    initGEOS(lwpgnotice, lwgeom_geos_error);

    g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
    if ( 0 == g1 )   /* exception thrown at construction */
    {
        HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
        PG_RETURN_NULL();
    }

    g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
    if ( 0 == g2 )   /* exception thrown */
    {
        HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
        GEOSGeom_destroy(g1);
        PG_RETURN_NULL();
    }

    retcode = GEOSHausdorffDistance(g1, g2, &result);
    GEOSGeom_destroy(g1);
    GEOSGeom_destroy(g2);

    if (retcode == 0)
    {
        HANDLE_GEOS_ERROR("GEOSHausdorffDistance");
        PG_RETURN_NULL(); /*never get here */
    }

    PG_FREE_IF_COPY(geom1, 0);
    PG_FREE_IF_COPY(geom2, 1);

    PG_RETURN_FLOAT8(result);
}

I don't have much experience with GEOS or interfacing it with PostGIS. The only thing I found in GEOS is that in that in geos::algorithm::distance::DiscreteHausdorffDistance, it is said that there is an double orientedDistance () member function. I am not sure if the documentation is even current.

Can some with the knowledge suggest how to adapt the above code to obtain the directed (or one-sided) Hausdorff distance from GEOS?

--EDIT--

The PostGIS info from SELECT postgis_full_version();' is as below (Ubuntu 16.04).

  POSTGIS="2.3.3 r15473" GEOS="3.5.1-CAPI-1.9.1 r4246" SFCGAL="1.2.2"
  PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released
  2015/09/16" LIBXML="2.9.3" LIBJSON="0.11.99" (core procs from "2.3.2
  r15302" need upgrade) TOPOLOGY (topology procs from "2.3.2 r15302"
  need upgrade) RASTER (raster procs from "2.3.2 r15302" need upgrade)
  (sfcgal procs from "2.3.2 r15302" need upgrade)


texts: Different FIPS Codes Versions for Counties

texts:
I am a SQL guy writing queries for supply chain consultants. I know little about FIPS.

I found in our database, at least three versions of FIPS code for counties. The number of FIPS are from 3100 to 3300.

This caused many troubles, if I inner join (merge) on the FIPS code, but the FIPS are not the same across the entire database.

How to manage the different versions of FIPS code?

Should I do something like referencing?

What if county A in 2000 split into county AA and county AAA in 2010?

Can you introduce how do you manage such issue?

"at least three versions of FIPS code" means there are at least 3 different total number of FIPS code (for county). The smallest is about 3140, the largest about 3250.

If the number of FIPS is not same for table A and table B, it will cause trouble when you join table A and table B.


texts:
I am trying to make an online map with GeoServer.

I am making the map now in HTML with the OpenLayers library. I want to use a different base map because the default map from OpenLayers is more than I need.

I only want to see the countries, for example, not the roads and things like that.

Is this possible? And if so, how is this possible then?


texts:
I have a large set of data points that essentially represent the coastline of the world. I don't know the exact resolution but I'd say it's in the neighbourhood of 5-10m.

Imagine a set of points which represent a jagged coastline.  At a 1m resolution, we see all the detail and need 50 data points to correctly represent the coast, however, at 10km resolution the coast looks like a straight line and we would only need 2 points.

I'm writing a mapping application that will require the data to scale from "full-planet" resolutions to "street" resolutions. I need some help reducing the resolution of my data points, and subsequently the amount of data for wider resolutions.

I've been reading about Vector-Tiling and think this will be the best solution (I'm doing something like this already), however, at the low resolutions, my dataset it still huge meaning that rending it takes a [relatively] long time, where as the high resolutions are fast because the 'effective' dataset it small (subset of the whole).

I'm trying to determine how to take my large dataset and reduce the scale of it so that my 'full-planet' view's tile data is a manageable size.  How do I create the tile sub-datasets from the larger complete dataset?

I'm happy to use a tool but I'd rather do it myself in an effort to learn
how it's done.


texts: Trouble intersecting a polyline and contours (Hack index / stream length gradient index)

texts:
Working on creating polylines for measurement using SL index (hack index). So far I have the main trunk of a stream, and contours at 10m. What I am trying to do is

a) get the elevation at each intersection of the stream (polyline) and a contour for the full length of the stream

b)get the distance along the stream between each contour intersection

So far I've been able use spatial join to get the contact between the stream and the contour but the distance is off (checked by using measure tool on the polyline of the full stream and summing them in the attribute table)

Any suggestions on how to make this work?

Thanks!


texts:
I'm trying to create some layers from feature classes located in the geodatabase. The MXD is created, symbology on feature classes defined within ArcMap, but for some reason my symbology is not carried over to layers. Do you know why?

import arcpy
import os

cwd = os.getcwd()
gdbLoc = cwd + "/../Zion.gdb"
lyrLoc = cwd + "/../symbology/"
mxdLoc = cwd + "/../Assignment3.mxd"

arcpy.env.overwriteOutput = True

mxd = arcpy.mapping.MapDocument(mxdLoc)
df = arcpy.mapping.ListDataFrames(mxd, 'Layers',)[0]

for fc in arcpy.mapping.ListLayers(mxd, "", df):

    outlayer = str(fc)+"_lyr"
    layerfile = os.path.join(lyrLoc, str(fc) + '.lyr')
    arcpy.MakeFeatureLayer_management(fc, outlayer)
    arcpy.SaveToLayerFile_management(outlayer,lyrLoc+str(fc))

    print '{0} created.'.format(fc)


texts: Setting QGIS user defined coordinate system (code problem)

texts:
I'm trying to set custom CRS in QGIS v 2.16.3. (Serbia Gauss-Kruger Zone 7). I've got parameters from .prj file that I use when working with GlobalMapper:

PROJCS["Transverse_Mercator",GEOGCS["GCS_Geographic Coordinate System",DATUM["D_SRBIJA",SPHEROID["Bessel_1841",6377397.155,299.1528131060786]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["scale_factor",0.9999],PARAMETER["central_meridian",21],PARAMETER["latitude_of_origin",0],PARAMETER["false_easting",7500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

Can someone help me and translate this code into a form that uses the QGIS? This should be a modification of MGI 1901/Balkans zone 7:

+proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +y_0=0 +ellps=bessel +towgs84=682,-203,480,0,0,0,0 +units=m +no_defs


texts:
I want to use arcpy to download a JSON file form "NoiseTube.net" and convert the file to a feature class. However, i'm finding it difficult  to complete this task, since the json file does not contain "geometry" and "type" keys. I want to be assisted to complete my task. Attached, is my script
import arcpy
import os
from urllib import request
import json
import tempfile

url = "http://www.noisetube.net/users/11463/tracks/88383.json"
temp_dir = tempfile.mkdtemp()
filename = os.path.join(temp_dir, "Noise.json")
response = request.urlretrieve(url, filename)
json_file = open(filename)
json_raw = json.load(json_file)

arcpy.management.CreateFileGDB(r"Z:\ARCGIS", "Noisetube.gdb")
arcpy.env.workspace = r"Z:\ARCGIS\Noisetube.gdb"

Example of the JSON:
{"dataset":[{"lat":null,"lng":null,"db":"40.136","date":"2019-12-08T20:19:02Z","id":64386058},{"lat":5.1178671771186535,"lng":-1.2908399999892466,"db":"37.85","date":"2019-12-08T20:19:03Z","id":64386059},{"lat":5.1178671771186535,"lng":-1.2908399999892466,"db":"36.986","date":"2019-12-08T20:19:04Z","id":64386060},{"lat":5.1178671771186535,"lng":-1.2908399999892466,"db":"37.483","date":"2019-12-08T20:19:05Z","id":64386061}]}


texts: Terrain Ruggedness Index (TRI) calculated in QGIS is approx 4 times that calculated in GDAL

texts:
When I calculate TRI using QGIS and GDAL tools for the same elevation DEM, I get results that are comparable between the two except that the absolute values in QGIS are approximately 4x that in GDAL. I am not modifying any defaults.
Can anyone shed light on what may be going on?
Which set of TRI values may I trust?


texts: MapServer dumps error: msPostGISLayerWhichShapes(): Error (ERROR: Operation on mixed SRID geometries)

texts:
Using the following WFS: http://dev-maps.isric.org/mapserv?map=/map/wosis_latest.map on QGIS 2.18 and then opening the QGIS query builder on layer ms:wosis_latest_profiles, and running the following query:
profile_id IN (69728,69916,69525,69903,69486,69607)
We end up with the following error on MapServer's side :
MapServer outputs a msPostGISLayerWhichShapes(): Error (ERROR:  Operation on mixed SRID geometries) error.
From the logs we have a rather extensive SQL query: https://pastebin.com/cQDXM0Hp
Basically MapServer sends a query to postgis that contains:
((st_intersects(geom,ST_GeomFromText(
'POLYGON ((-184.2090000000000032
-184.2090000000000032 -54.8890999999999991))')

And here is the problem, for the query to be successful we need for ST_GeomFromText to indicate the SRID=4326, otherwise we have the error  Operation on mixed SRID geometries and no reply from MapServer
The WFS service runs mapserver7.6/gdal2.4/proj4 and postgis ("POSTGIS="2.5.4" [EXTENSION] PGSQL="110" GEOS="3.8.1-CAPI-1.13.3" PROJ="Rel. 7.0.1, May 1st, 2020" GDAL="GDAL 3.0.4)
The mapfile has the  indication of layer projection and simple data connection:
 LAYER
        NAME "wosis_latest_profiles"
    PROJECTION
        'init=epsg:4326'
    TYPE POINT
    CONNECTIONTYPE POSTGIS
    CONNECTION "...."
    PROCESSING "CLOSE_CONNECTION=DEFER"
    DATA "geom from web_services.wosis_latest_profiles"

Is this a problem related to the PostGIS version / GDAL combinations, some miss configuration on MapServer, bug on QGIS? Any tips on how to debug it?


texts:
I'm doing a research project that involves using demographic data collected from both the national Census and the American Community Surveys.

I am interested in about 40 zip codes within NYC.

The tool I was using to query Census and ACS only can pull up the most recent year, but I need to go a few years back.

What tool should I be using?

For example, I'm looking for the estimated population between ages 35 to 44 in each of 40 zip codes in 2013 as well as 2011.


texts: Is there any geocoding for better international support?

texts: I need international support database with Chinese and English, do you know any data provide both? (So far, I found only Google Map and MaxMind is okay;  Open Street Map just poorly supported data; any better solution?)


texts: Creating a source and a layer from vector tile features in Openlayers

texts:
So I retrieve all my features from a vectorTile layer (PBF from Geoserver), and now I want to recreate a layer by filtering on those features and getting only those that satisfy a request.
I successfully stored my filtered features in an array, but when I try to add them to a new vector source, I get the error: typeError e is null.

See my code below:

var filteredFeatures= [];

for (var i = 0; i < data1.length; i++) {
     for (var j = 0; j < features.length; j++) {

         if (features[j].get(data1[i].type) == data1[i].valeur) {

             filteredFeatures.push(features[j]);
         }

     }

     var source = new ol.source.Vector({
         features: filteredFeatures,
         projection: proj2154,
     });
      map.addLayer(new ol.layer.Vector({
          source:source,
           style:styles.red
      }));

     filteredFeatures = [];

 }


texts: Unknown coordinate system for the Netherlands

texts: I have a table with x and y coordinates. It is mentioned that coordinates in "RD Bessel" system. I couldn't find it in ArcGIS and QGIS. Also no info on google?


texts: Creating PostGIS Geodatabase features in QGIS defaults to ST_Geometry

texts:
I have a PostGIS database that uses the PG_Geometry stored in a shape field. This same database is a registered ArcGIS geodatabase. I am able to edit and create features in ArcGIS. I can connect, view and edit attribute data in QGIS, but when I try to create new features I get the error shown below. Ironically it seems that QGIS attempts to use ST_Geometry when writing the data even though the data are stored as PG_Geometry. I am looking for the setting in QGIS that will allow me to write the geometry data using PG_Geometry. Any help is much appreciated.

Could not commit changes to layer drp_occurrence

Errors: ERROR: 1 feature(s) not added.

Provider errors:
      PostGIS error while adding features: ERROR:  column "shape" is of type geometry but expression is of type st_geometry
    LINE 1: ...atex","dgupdatey","catalognumbernumeric") VALUES (st_geomfro...
                                                                 ^
    HINT:  You will need to rewrite or cast the expression.


texts: I use JavaScript to work with the program OSM2PO. I want to solve the problem of TSP. I need to get the order Vertex Ids and array waypoints in one query. Is it possible in the format Json? It would also be good to get the distance and time.


texts:
When developing an ArcMap extension, is there an automatic way to handle versioning?

In the snippet below, the Version tag can be updated manually, but is there a way to auto-increment it or let the build process version it?

<ESRI.Configuration xmlns="http://schemas.esri.com/Desktop/AddIns" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">  <Name>MyExtension</Name>  <AddInID>{...}</AddInID>  <Description>My Extension</Description>  <Version>1.0</Version>  <Image>Images\MyExtension.png</Image>  <Date>13/03/2013</Date>  <Targets>    <Target name="Desktop" version="10.0" />  </Targets>  <AddIn language="CLR4.0" library="MyExtension.dll" namespace="MyExtension">  </AddIn></ESRI.Configuration>


texts: How can I get real time log in .bat file at QGIS & PlugIn with Python PyQT4

texts:
I try to get real time log in batch file with python, then show it to QGIS Plugin real time. ( .bat file have run -jar file )

I success logging but not real time log in QGIS because QGIS is stopped by .bat file ( it is not completed )

try:
    process = subprocess.Popen(
        'C:/val/start.bat',
        shell=False,
        stdout=subprocess.PIPE,
        #close_fds=True,
        #creationflags=DETACHED_PROCESS
    )

    while True:
        line = process.stdout.readline().decode('utf-8')
        self.dockwidget.plainTextEdit.appendPlainText(line)
        QgsMessageLog.logMessage(line, tag="Validating", level=QgsMessageLog.INFO)
        if not line : break

    self.iface.messageBar().pushMessage("Validating", "Validating is completed")

except:
    self.dockwidget.plainTextEdit.appendPlainText("##############Error##############")

I want not pause QGIS on running .bat file and logging at QGIS Plugin by real time.

How can I do that??
Is it possible??


texts: Boolean test for file geodatabase having already been compacted?

texts:
We receive many 10.3 file geodatabases (fGDB) from a vendor. One contractual requirement is that the fGDB be compacted. Understandably we quality assure everything that comes in our door. Currently I have to record the file size of the fGDB and then run a compact, and look at the file size again. I can think of several ways of doing this with arcpy modules. Is anyone aware of a way to test for compaction without compacting the fGDB?. A true false return in ArcPy return would be handy.

I've not found anything in the forums and I'm inclined to think I might have to lodge an ArcGIS Idea.


texts: How to create a linestring from a geopandas data frame?


texts: ArcGIS XML Import Failed - how to find out why it failed?

texts:
I am creating an ESRi XML Workspace document in some of my own C# code. The XML doc that is produced validates against the GdbExchange.xsd schema provided by ESRI. However, when I use ArcCatalog to import the XML document into a newly created Personal Geodatabase (.mdb) I get the following error:

Import xml data failed
FDO error: -2147186428 [XML parse error: Unspecified error
 [Line: 415, column: 23].]

The line and column indicated above are:

    ...
    </DatasetDefinitions>
    <Metadata xsi:type="esri:XmlPropertySet">
        <XmlDoc />
    </Metadata>
</WorkspaceDefinition> <!-- THIS LINE -->
<WorkspaceData xsi:type="esri:WorkspaceData">
    <DatasetData xsi:type="esri:TableData">
        <DatasetName>PointInspections</DatasetName>
        ...

(Column 23 is after the closing angle bracket, ">", an the line marked "< !-- THIS LINE -->")

Given that the XML document is valid, how can I get more useful information from ArcGIS?

(There is a somewhat similar question asked here: Importing XML file in ArcGIS failed, but that didn't have a resolution that worked for me...)


texts: Calling a plugin in QGIS Graphical Modeler

texts:
I am familiar with Arcmap modelbuilder but struggling with QGIS equivalent. I have installed plugin 'UMEP' and need to run 'sky view factor' module on series of rasters.

Inside the Graphical modeler, I can only see modules which were by installed by default and not additional plugin.

Is there anyway which to drag and drop plugin in Graphical Modeler, like you would do with any Toolbox in arcmap?


texts: Creating Heatmap Based on Distance From Points

texts:
currently with postGIS I have created a table of point geometries at different coordinates representing cell phone towers. I would like to create a heatmap which would show the strength of the cell phone signal across the region where towers are located. So for example, create a grid and at each point in the grid first find all the cell towers within a certain radius of my grid point, then calculate the distance to each of the tower, then using this distance find a signal power value from a function based on distance, and finally plot a color value representing signal power based on combining all the distances to the points.

Is this possible in either postGIS or geoserver?


texts: I have a PostGIS table with linestring geometries and point (A) on it and I need to detrimine points (M) which located M meters from A. Distance is measured along linestrings. Any thoughts?


texts: Remove certain substrings from name column

texts:
I have 2 kind of Strings in column "NAME":

Berlin
Reg.-Bez. Hannover

I need a new colum "NAME2" changing the cities to:

Berlin
Hannover

I have this so far:

if("NAME" LIKE 'Reg.-%',xxxx,"NAME")

Can you help me with "xxxx"?

I tried

xxxx= "right( "field", length("field" ) - 10)"

but that's not working.

Or can i just use something like

trim('Reg.-Bez. ')


texts: Clipping Single Raster with Multiple Polygons in ArcGIS ModelBuilder?

texts:
I am quite new to using iterators, and ModelBuilder in general. I have a .tiff format raster (not a feature class) and I would like to clip it with 18 polygon shapefiles which are spatially present within the raster extent.

I looked at some other model formulas, such as Batch Clip in ModelBuilder of ArcGIS for Desktop?

But, I just can't seem to get it working in ArcGIS 10.2.

Any ideas how I can set up a model which would automate this process and export 18 unique clipped rasters?

I have edited below and added the final correct model.


texts:
Can somebody help me with gdalwarp issue? I am having an input raster and would like to clip it with polygons. After running the below code, while overlaying output raster on input raster, I see there is a slight shift in pixels of output raster.

Can somebody help me how to match the output and input raster so that the pixel borders matches exactly?

import ogr
import subprocess
inraster = 'C:/Users/claudio/workspace/test/inraster.tif'
print(inraster)
inshape = 'C:/Users/claudio/workspace/test/country.shp'
ds = ogr.Open(inshape)
lyr = ds.GetLayer(0)
lyr.ResetReading()
ft = lyr.GetNextFeature()
while ft:
country_name = ft.GetFieldAsString('admin')
outraster = inraster.replace('.tif', '_%s.tif' % country_name.replace(' ', '_'))
subprocess.call(['gdalwarp', inraster, outraster, '-cutline', inshape, '-crop_to_cutline', '-cwhere', "'admin'='%s'" % country_name])
ft = lyr.GetNextFeature()
ds = None


texts: Sending compressed folder with shapefiles and QGIS?

texts:
When I compress my QGIS file with the folder containing all my shapefiles, how can I maintain the paths?

For example, I made a map on QGIS using a large amount of shapefiles. Once done, I save my QGIS file to
C:\Users\Workcomputer\Desktop\New_folder

From there I compress the shapefiles and the QGIS file located in \New_folder and send the newly made .zip file to myself by email or dropbox.

However, when I download and uncompress the .zip on a different computer, the QGIS file is empty, and I have to re-path all the shapefiles, because they are now located in

C:\Users\Homecomputer\Desktop\New_folder

instead of

C:\Users\Workcomputer\Desktop\New_folder

Is there anyway to short-cut having to do this?


texts: Checking if attribute name exists using QGIS Field Calculator

texts:
I'm trying to write and reclassification expresion, within the field calculator, that works in two different layers (vectorial) with two different names for ther same thing. The issue here is that in my approach, as one of the attributes names doesn't exist in the layer, the expression doesn't work.

case
when ("layer1colname" = 10 or "layer2colname" = 10) then "CAT10"
end

I am  thinking about two different approaches:

Including a piece of code that first finds if the attribute exist.
Use other reclassification tool (for vectorial data) out of the field calculator.


texts: When to choose MapScript over a Map file or vice versa?

texts: Are there things that can be done with MapScript and not Map file or vice versa?  What are the scenarios to use each one?


texts: How to extract data points based on date attributes

texts:
I have around 1.5 million data points of fire incidents during a period of 10 years. I want to extract each month into a new layer so then I will basically have 120 layers in total.

Is there a way to do this without doing it manually one by one i.e. selecting points of every month and extracting it 120 times?

EDIT:
Software: ArcGIS 10.4.1
Data: shapefile


texts: Freezing legend on map composer so that it stops updating when changing symbols on map?

texts:
This seems like a counter intuitive question, but is there a way to freeze the legend so that it stops updating itself in the print composer?

My auto-update box is unticked for the legend, but the legend symbols keep updating with every map change/refresh. I want the symbol sizes to stay the same regardless of what I change in my main map.

I've tried to save the legend as an image and re-importing the image as the legend, but the image is not scalable and therefore looks very pixelated in comparison to the rest of the map.

I am currently using Qgis 2.4


texts:
I have this coordinate reference system defenition:
GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],VERTCS["NAVD_1988_Foot_US",VDATUM["North_American_Vertical_Datum_1988"],PARAMETER["Vertical_Shift",0.0],PARAMETER["Direction",1.0],UNIT["Foot_US",0.3048006096012192]

I run this test using osr python module to check if projection definition is valid.
import osr
srs = osr.SpatialReference()
srs.ImportFromESRI(['GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],VERTCS["NAVD_1988_Foot_US",VDATUM["North_American_Vertical_Datum_1988"],PARAMETER["Vertical_Shift",0.0],PARAMETER["Direction",1.0],UNIT["Foot_US",0.3048006096012192]]'])
srs.Validate()

This outputs 5, which is OGRError for corrupt data.
When I use this same string to validate crs definition in QGIS it is saying the proj definition is valid.

How can a spatial reference system be corrupt per OSR while valid per QGIS?


texts: Summarize forest change by geometry, year

texts:
I am trying to calculate the forest change from the Hansen Global Forest Change dataset on Google Earth Engine.

My end goal is a columnar output of forest change by year by county.

I was able to calculate the total net change from 2000 to 2018 using this tutorial, but I cannot figure out how to break it out by year.

How do I summarize/reduce by a geometry AND a year? Do I have to loop through years?


texts: Update Layer Symbology based on 1 field with changing data

texts:
I am running ArcGIS 10.2.2 Version 10.2.2.3552, using ArcMap.

I have quite a big task here and before i attempt it i want to know if it is even possible.  I've looked up some similar questions on here and have some bases, but I need more of a directed starting point.  I have several mxds (100 or so) and each month the data changes.  I have 15 colors and 1 grey color for "other" in my symbology based on 1 field.  Each month the main colors could change though, based on if they had more data then the previous month.

I export these maps based on data driven pages with python.  Is there a way to update symbology based on new statistics for my data driven pages? So the colors are mainly the same but a few rotating colors every month and the old ones going grey, new ones taking their place.  I get the statistics through the "Summary" tool, but i can't find an equal python function for this.

The thing I've looked into the most is representation.  Though I am not sure of the limitations of these and if these even will work with what I am trying to do?


texts:
I kept a listing of virtual globes in 2005 and 2006, but it's now woefully out of date http://geography2.blogspot.com/2005/11/listing-of-virtual-globes.html. Today, the main four seem to be: Google Earth (proprietary and free), Microsoft Bing Maps 3D (proprietary and free), ESRI ArcGIS Explorer (proprietary and free), and NASA World Wind (free and open source). Many, if not most, are being integrated into the web browser.

What is out there? Any special or particularly innovative capabilities?


texts: QGIS Importing several CSV-files using Python

texts:
I am trying to import several CSV-files using Python. However, my code does not work. Furthermore, the coordinates are in DMS-format and I am unsure on how to implement that to my code.
import os.path, glob
layers=[]
for file in glob.glob('/Users/NAME/Desktop/gis/weights/*'):
    uri="file:///"+file+"?type=csv&xField=LON&yField=LAT&spatialIndex=no&subsetIndex=no&watchFile=no&crs=epsg:4326"
    vlayer = QgsVectorLayer(uri, os.path.basename(file), "delimitedtext")
    vlayer.addAttributeAlias(0,'x')
    vlayer.addAttributeAlias(1,'y')
    layers.append(vlayer)
QgsMapLayerRegistry.instance().addMapLayers(layers)

When running my code I get the following error-messages:
Traceback (most recent call last):
File "/Applications/QGIS3.10.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/code.py", line 90, in runcode exec(code, self.locals) File "<input>", line 1, in <module>
File "<string>", line 10, in <module> NameError: name 'QgsMapLayerRegistry' is not defined

Below is a sample of the data in the CSV-file.


texts: Defining uphill side of roads using ArcGIS Desktop?

texts:
I have a big raster of a mountainous region, and I need to extract the uphill side of roads, from the road to the peaks of the mountains.

For a road, I know the uphill side changes to become the downhill side when you make a hairpin turn. For this reason, I believe the processing should start from the lowest road and continue up into the mountains (the highest points).

Is there any possibility of doing this?

I don't have a good knowledge of the Python language.


texts: Deleting close points from FeatureCollection using buffer in Google Earth Engine?

texts:
I want to generate a number of random points over a given area, e.g. all land masses of the world. For this I have coded

var fcLandMasses = ee.FeatureCollection('ft:1l9VK5FJ4_hiH2gHjr0MtxwB-CrNEqOE5EBlW6Q');
var fcRandomPoints = ee.FeatureCollection.randomPoints(fcLandMasses,numberOfPoints,seed);

where you need to set numberOfPoints to the number of points you want to generate and seed to an arbitrary non-negative integer to ensure repeatability.

How can I clean up the FeatureCollection fcRandomPoints such that points located too close to each other (i.e. within a certain buffer) are removed?

I think there are two general choices:

Keeping one (e.g. the first-appearing) point of points within that buffer.
Keeping the average of all points within that buffer.

I'd like to follow option 2, which will probably include a mapping over fcRandomPoints.

What is the most efficient way to implement this?


texts: How to update popup after dragging the marker?

texts:
I have a leaflet map where i am using Control Geocoder to get some locations on search,then i add the coordinates on 2 input fields,also i added a popup for the search result,the problem is that when the marker is moved the popup doesn't update with the new place,how can this be done?

Mmap = L.map('Modalmap').setView([8.7144, 125.7481],8);
L.tileLayer( 'http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    attribution: '&copy; ',
    subdomains: ['a','b','c']
}).addTo(Mmap);
var geocoder = L.Control.geocoder({
    defaultMarkGeocode: false,
    draggable:true
})
.on('markgeocode', function(e) {
    var box = e.geocode.center;
    var name=e.geocode.name;
    document.getElementById("Latitude").value = box.lat;
    document.getElementById("Longitude").value = box.lng;
    MarkLayer=L.marker([box.lat,box.lng],{draggable:true}).addTo(Mmap).
    on('dragend', onDragEnd).
    bindPopup(e.geocode.name).
    openPopup();
    displayLatLng(box);

    group = new L.featureGroup([MarkLayer]);

    Mmap.fitBounds(group.getBounds());
}).addTo(Mmap);

function onDragEnd(event)
{
    var latlng = event.target.getLatLng();

    displayLatLng(latlng);
}

function displayLatLng(latlng) {
    document.getElementById("Latitude").value = latlng.lat;
    document.getElementById("Longitude").value = latlng.lng;
}


texts: Accessing slope and aspect data using Leaflet?

texts:
I am building an application in leaflet, and I need to be able to access slope and aspect as a function of latitude and longitude.  Currently I am able to get elevation data from the mapbox elevation dataset here, which can be processed pretty simply to get elevation as a function of latitude and longitude.

Does anyone know of a similar dataset, accesible in leaflet, that gives slope and aspect?

I know that ArcGIS online has a terrain dataset which contains these things.  There is also a leaflet esri plugin that allows you to add the terrain layer in with L.esri.basemapLayer('Terrain').  But I don't see a way to actually extract data from that, just to view it within leaflet.

There is this article which praises leaflet's creator for his ability to calculate aspect and slope in realtime as the user pans or zooms, based on the mapbox elevation dataset I mentioned above.  But frankly I don't understand the code that he's quoting, and this seems both inaccurate at low zoom levels, and computationally expensive for high resolution maps.  I thought of potentially creating my own function which takes a point, gets 6 or 8 points around it, get the elevation of those points, and then calculates slope and aspect from that.  But that seems inefficient and computationally expensive, especially when I'm going to need to get slope and aspect for hundreds of points many times per second.

I have been casually searching around for months for how to get slope and aspect data within leaflet.  If I have to learn a whole new framework to import into leaflet, I'm open to that.  I'm just not sure where to start.  Forgive me if this question is in the wrong place or has been answered elsewhere.


texts: Making coordinates match between Google Maps and matplotlib.basemap?

texts:
I'm trying to plot some data with matplotlib and basemap, but I have trouble making the coordinates match between Google Maps and basemap. For example, I click at a distinctive feature on the map in Google Maps, and read the coordinates (this is on the north east coast of Greenland:

I then try to plot those same coordinates on a map with matplotlib and basemap:

The two points seem to differ by about 10 km. Are there any parameters I can pass to basemap, or any conversions between ellipsoiods or anything I can do to make these points match up?

The python code I used was this:

#!/usr/bin/env python

import matplotlib
matplotlib.use('agg')
from mpl_toolkits.basemap import Basemap
import mpl_toolkits.basemap.pyproj as pyproj
import matplotlib.pyplot as plt

# These are the coordinates I read out from Google Maps
lat, lon = 80.615439, -15.474109

m = Basemap(width=200000, height=150000, resolution='i', projection='stere', lat_ts=80, lat_0=lat, lon_0=lon)

m.fillcontinents(color='#aaaaaa', lake_color='#cccccc')
m.drawcountries(linewidth=0.2)
m.drawmapboundary(fill_color='#cccccc')

X, Y = m(lon, lat)

plt.scatter(X, Y, s = 20, color = 'red')
plt.savefig('map.png')

I tried this:

# These are the coordinates I read out from Google Maps
lat, lon = 80.615439, -15.474109

gmaps = pyproj.Proj('+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs')
wgs84 = pyproj.Proj(init='epsg:4326')

print(lat, lon)
lon, lat = pyproj.transform(gmaps, wgs84, lon, lat)
print(lat, lon)

This converts from 80.615439 -15.474109 to 0.000724 -0.000139.

What I'm looking for is a way to convert between two different lat-lon systems.


texts: Could a QGIS plug-in conceivably save and restore the main-window geometry of any project?

texts:
With respect to this:
Restoring the exact original situation of a spatial bookmark?
and this:
How to Resize the Qgis Main Window from python
Q. Could a plug-in conceivably save the current state of the main QGIS window so that it could be restored later prior to navigating to a saved spatial bookmark?
It would make it much easier to add and print new layers in old projects which would be guaranteed to line up with previously printed layers. (I know that layers would align in QGIS, but I refer to the printed-out images registering correctly.)
The question could be replaced with:
Q. If I have saved a spatial bookmark for a project, how can I add a layer to that project that will be guaranteed to align with layers printed earlier?
(I had the sense from the first reference that it was not possible to completely, automatically, reconstruct the earlier window geometry at a later point. Indeed, I've now tried it and found the edge-positioning of the window to be too coarse to get perfect reconstruction.)


texts: ArcGIS JS API - Can't draw results of Query Task

texts:
I have an ArcGIS JavaScript map that I'm trying to use with a Query Task. The query runs fine on the server side and the proper geometry is returned. Then, I draw the feature as follows:

var features = featureSet.features;

for (var i = 0; i < features.length; i++) {
    var feature = features[i];
    feature.setSymbol(lbSymbol);
    feature.setInfoTemplate(lbInfoTemplate);
    map.graphics.add(feature);
}

This was working for a while but it randomly broke and I can't figure out why. No errors -- the map just doesn't draw anything. If I listen for the graphics-add event on the graphics layer it fires and all the event data looks fine (including spatial reference, geometry, etc.) Can anyone think of why this wouldn't work or how else to troubleshoot?

For what it's worth - this code to zoom to the extent of the feature is also broken:

var extent = esri.graphicsExtent(features);
if (extent) {
    extent = extent.expand(3.0);
    map.setExtent(extent);
}

The extent of esri.graphicsExtent(features) matches the extent of the feature. Not sure why the map doesn't zoom on setExtent.

UPDATE: I noticed that graphics layer is referenced to WGS84, while the map and graphic feature are State Plane. I'm just using the built-in map.graphics layer -- I assumed this would match the spatial reference of the map automatically. Still puzzled.


texts: How to tune the performance of Nominatim Reverse Geocoding Application

texts: Running the query directly to the PostgreSQL takes 3-8ms. The PostgreSQL can handle 500 concurrent connections. The Apache server is running on 8 core CPUs. Analyzing the Apache access logs: There are 14-21 request served per sec. Requesting via the browser, it takes about 300ms TTFB waiting time, 3ms content download, 1ms initial connection (ping). I would really want to know the time spent to convert the request to a query and also the time taken to process the data from the PSQL to the response. What should I do to reduce the waiting time (TTFB)?


texts: WMS layer above select-controlled vector layer in OpenLayers 2?

texts:
Is it possible to display an image layer (WMS or Image) above a vector layer with activated SelectFeature control in OpenLayers?

The image layer (radar returns) contains only sparse, semi-transparent regions and doesn't significantly obscure the underlying vector layer, whose features I would like to be able to select. It seems reasonable to expect this behaviour from the map.

However, as soon as the OpenLayers.Control.SelectFeature control for the vector layer is activated, the vector layer is moved to the top of the map, placing it above the image layer.  Calls to map.setLayerIndex() have no effect.

Among many attempts to fix this, I was able to prevent it from occurring by omitting the moveLayerToTop() function in the OpenLayers code, as suggested here:

https://stackoverflow.com/questions/4728852/forcing-an-openlayers-markers-layer-to-draw-on-top-and-having-selectable-layers

Although the image layer now appears above the vector layer, the OpenLayers.Control.SelectFeature no longer functions whenever the image layer is set to visible.

The image layer appears to be blocking mouse events to the control, although the map itself responds normally.  I tried setting the image layer's events.fallThrough to true but it had no effect.  I also (in desperation)
set the layer's pointer-events to none; again, to no effect.

I saw this thread:

OpenLayers how to maintain layer ordering while using SelectFeature control

and was led to the featureclick.js extension:

https://github.com/openlayers/openlayers/pull/174

I'll likely give that a try next, but am wondering if I've missed something obvious in standard OpenLayers.  I haven't found any other solutions posted on the web.


texts: UpdateGeometrySRID and ST_Transform are the same thing?

texts:
UpdateGeometrySRID actually transforms the coordinate system of my data?
For in the PostGIS documentation it appears the trim ST_Transform (Returns a new geometry with its coordinates transformed to the SRID referenced by the integer parameter.).
What is the difference between these functions?

Another thing the sql below really redesigns same data? So no need to use UpdateGeometrySRID

ALTER TABLE lotacao ALTER COLUMN geom TYPE geometry (MULTILINESTRING, 101020) USING ST_Transform (ST_SetSRID(geom, 101010), 101020);


texts: 'GIS' grade mobile mapping app that isnt ArcPad?

texts:
I am a GIS analyst/geologist with a mining company based in Canada. However most of my fieldwork is in remote places with minimal access to the internet. I am after a reliable alternative to ArcPad and ESRI running on either ipad or android tablets. I have already tried several apps like Mapitfast and cartomobile but Ive been unimpressed with either their UI and/or inability to share/access data off the device.

Some tools I need are:

Creating rock outcrop maps
Editing/importing geology polygons
GPS tracking that can be exported for later use
Photo-tagging would be a plus
Bing or Google as a background image

Thanks in advance.

-David-


texts: Is it possible to import a animated GIF into QGIS? And is it also possible to georeference it, to have it on the right place?


texts: Select point(s) from given coordinate and move to new coordinate using arcpy

texts:
Goal: Move address points to their correct location.

Plan: I have a layer of AddressPoints of which some need moved to a new location.  I want to take a line layer (called AddressPoints_AdjustLocation, shown as a purple dashed line below) and draw a line from the address points current location (which is usually close to the road) to the desired location (which is the front door of the building)

Code I have so far:

    arcpy.MakeFeatureLayer_management("../pathToAddressPoints", "lyr_AddressPoints")
    arcpy.MakeFeatureLayer_management("../pathToAddressPoints_AdjustLocation", "lyr_AddressPoints_AdjustLocation")

    # Create a cursor that holds all the selected AdjustLocation features
    rows = arcpy.SearchCursor("lyr_AddressPoints_AdjustLocation")
    if rows:
        for row in rows:
            feat = row.getValue("SHAPE")
            vertex_list = []
            for vertex in feat.getPart(0):
                if vertex:
                    vertex_list.append([vertex.X, vertex.Y])

            vertexStart = vertex_list[0]
            # TODO: Select all address points in "lyr_AddressPoints" that intersect the vertexStart location

            vertexEnd = vertex_list[-1]
            # TODO: For each address point selected above, move it to the vertexEnd location

The two items I need help with are how to select address points that Intersect the vertexStart position, and then moving those address points to the vertexEnd location (using UpdateCursor I assume).


texts: Adding Encryption options to GeoServer SQL Server connection

texts: I have set up an HTTPS connection between GeoServer and an Azure database using the SQL Server extension, however I would like to add encryption options to the connection string.  Is this possible in GeoServer 2.16 (currently running on Tomcat 9)?  Would I have to modify the actual extension itself?


texts: QGIS Server no fonts getting rendered on Ubuntu 18.04

texts:
I've got my QGIS Server instance up and running and features are being rendered correctly via WMS. I can't get any labels to show up, however. I'm sure my problem is related to this post, but i'm not getting the blue boxes showing the missing fonts at all.

-I'm running Ubuntu 18.04 and installed it based on this great tutorial

-I can't see anything in the qgis server error log because it is not getting written to. This is another thing I can't figure out but it has the same privs as the apache2 error log:

-i've run:

export QT_QPA_FONTDIR=/usr/share/fonts/truetype/dejavu

-I've added the SetEnv to my apache qgis-server.conf:

-I'm pretty sure that i'm referencing the correct font on the system and in the .qgs file. $fc-list shows that DejaVu Sans is loaded and the style name is "Book"

And the reference in the .qgs file is:

So, i'm totally stuck. Any ideas on where to hunt next or why this font isn't working? I'm also not sure if i'm missing the QT package all together? It doesn't appear to be in aptitude.


texts:
I would like to display an ImageWMS layer using OpenLayers v6.3.1, OL-Cesium v2.9 and Cesium v1.70 in 3D mode. In 2D mode it is perfect, however in 3D mode, I don't see the states ImageWMS layer, the map show only the OSM basemap. There are no errors on the browser's console and on the Network tab there are no WMS requests.
Any suggestions? How can I solve this issue?
This is my code below:
<head>
    <link rel="stylesheet"
        href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.3.1/css/ol.css">
    <script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.3.1/build/ol.js"></script>
    <link href="https://cesium.com/downloads/cesiumjs/releases/1.70/Build/Cesium/Widgets/widgets.css" rel="stylesheet">
    <script src="Cesium.js"></script>
    <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/ol-cesium@2.9.0/css/olcs.css">
    <script src="olcesium.js"></script>
    <style>
        html,
        body,
        #olCesMap {
            height: 100%;
            width: 100%;
            padding: 0;
            margin: 0;
        }

        #btnSwitchMap {
            position: absolute;
            top: 10px;
            right: 50px;
            z-index: 2;
            padding: 4px;
        }
    </style>
</head>
<body>
    <div id="olCesMap">
        <button id="btnSwitchMap" onclick="cesViewer.setEnabled(!cesViewer.getEnabled());">Switch 2D/3D</button>
    </div>
    <script defer>
        var controls = ol.control.defaults({
            rotate: true,
            rotateOptions: {
                tipLabel: "Reset rotation. \nUse Alt+Shift+Drag to rotate the map."
            }
        });
        var interactions = ol.interaction.defaults({ altShiftDragRotate: true, pinchRotate: true });

        // Declaring Zoom level to further view
        var zoom = 5;
        // Declaring Map center
        var center = ol.proj.fromLonLat([-96.41778916767144, 39.90201978025539]);

        const basemapTileLayer = new ol.layer.Tile({
            title: 'BaseMaps',
            displayInLayerSwitcher: false,
            source: new ol.source.OSM({
                'url': "https://tile.openstreetmap.org/{z}/{x}/{y}.png"
            })
        });
        basemapTileLayer.set('name', 'tileLayer');

      const imageWMSSource = new ol.source.ImageWMS({
        url: 'https://ahocevar.com/geoserver/wms',
        params: {'LAYERS': 'topp:states'},
        ratio: 1});

        //Map Creation
        var map = new ol.Map({
            controls: controls,
            interactions: interactions,
            target: 'olCesMap',
            layers: [basemapTileLayer,new ol.layer.Image({
      extent: [-13884991, 2870341, -7455066, 6338219],
      source: imageWMSSource
    })],
            view: new ol.View({
                center: center,
                zoom: zoom
            })
        });
        //Cesium Viewer constructed based on Ol.Map
        var cesViewer = new olcs.OLCesium({
            map: map,
            sceneOptions: {
                mapProjection: new Cesium.WebMercatorProjection()
            }
        });
        //Cesium Viewer scene
        var scene = cesViewer.getCesiumScene();
        //Cesium Viewer enable
        cesViewer.setEnabled(false);

        //Cesium Viewer depthtest
        scene.globe.depthTestAgainstTerrain = true;
    </script>
    </body>

UPDATE: I tried with Ol-Cesium v2.10, but I got these errors:

Uncaught SyntaxError: import declarations may only appear at top level
of a module olcesium.js:4 Uncaught ReferenceError: olcs is not defined


texts:
I have code written in ArcMap which uses the package ArcPy which I cannot use in a stand-alone Python script, but I want to do it with GDAL without ArcMap. Basically it goes through the lines in the shapefile and create points every 20 meters and stores all these points in a new shapefile.

I have tried using the interpolate method in fiona and am just very lost.

Could someone help?


texts: My point cloud orientation is very wrong

texts:
I'm trying to load a cloud point using Cesium and Entwine.

This is my code until now:

First I convert a LAS file to EPT using this tutorial but need to convert projection to EPSG:4978.

entwine build -i https://data.entwine.io/red-rocks.laz -o ~/entwine/red-rocks-ept -r EPSG:4978

Now I have an EPT folder. Lets convert it to Cesium 3D tiles:

entwine convert -i ~/entwine/red-rocks-ept -o ~/entwine/red-rocks

Now I have a folder with Cesium files.

Lets serve them. I'll use docker but any http server will do it well. Don't forget to allow CORS:

docker run -it -v ~/entwine/cesium:/var/www -p 8080:8080 \
    connormanning/http-server

Now in Cesium:

var tileset = scene.primitives.add(new Cesium.Cesium3DTileset({
     url : 'http://YOUR_SERVER:8080/red-rocks/tileset.json',
}));

tileset.readyPromise.then(function(tileset) {
    camera.viewBoundingSphere(tileset.boundingSphere, new Cesium.HeadingPitchRange(0, -0.5, 0));
    camera.lookAtTransform(Cesium.Matrix4.IDENTITY);
});

My cloud is floating or lost somewhere in space:

Anything I do with mouse (change the camera view) the cloud disapear forever.

Already tried to remove the readyPromise event.

Can't find the ENTWINE tag.


texts: Showing only some layers in QGIS legend?

texts:
My QGIS 2.8.2 standalone app provides several decorating layers (e.g. cartographic grid, background layers, reference points) which should be not changed by the user (neither editing nor hiding), as well as layers, the user can modify and query.

The app has a legend widget which shows all layers, since all layers are registered with map registry. The widget is populated with the layer tree as described in How to add a legend to a canvas in a standalone PyQGIS application?.

Instead of exposing all layers in the tree view, I want to restrict the tree view to the layer group, whose children represent the user mutable layers. Thus the user can change visibility and ordering of only these layers, and must keep the decoration as it is.

What I tried so far is to initiate QgsLayerTreeMapCanvasBridge with this group node (instead of layerTreeRoot), pass the group node to QgsLayerTreeModel(), or
changed the root group of the model with QgsLayerTreeModel.setRootGroup(). In each case Python crashes, with an error message like this:

QObject::connect: Cannot connect (null)::willAddChildren( QgsLayerTreeNode*, int, int ) to QgsLayerTreeModel::nodeWillAddChildren( QgsLayerTreeNode*, int, int )

Is it possible to expose only some group oder layer nodes in the treeview, and how can I implement this?


texts: How can I install a python 3 module for its use in QGIS3 in Windows?

texts:
I've been struggling to install BeautifulSoup4 for a Python script I want to share with Windows users, but I just can't get it installed for its use in QGIS,
I've tried

pip install BeautifulSoup4

and using the ez_setup.py

in the OSGeo4W, but python --version shipped with osgeo4w is 2.7.14 (latest osgeo4w release), while QGIS uses python 3, so I can't even use pip3 because osgeo4w shell is 2.7.14;

Script works good on linux.


texts:
I'm trying to identify the mountain roads in a route, and I've thought that a good approach could be get the paths with a very closed curve (like a U turn), but now I can not find any documentation about.
Do you know any idea, approx or some documentation or article where this topic has been covered?
The main idea is only know if a route has curve like the image below, or not. For this in my app I'm using PostGIS and Python and data (rote path) comes from OSM.


texts: Converting JP2 to GeoTIFF QGIS 3.16 Win 64

texts:
I have hundreds of USGS NAIP imagery GeoTIFFs that are in JP2 format, EPSG:3857.
I am using QGIS 3.16 to convert them to GeoTIFFs ("Rendered Images") as EPSG:4326, but I can only process one image at a time from within the GUI which is very, very tedious...
I can't use the GDAL_TRANSLATE command line option as when I originally built GDAL quite some time ago, I did not foresee using JP2 imagery. I don't have the time right now to rebuild GDAL (and a little afraid I might break something and create a work stoppage...), so I'm using GeoJasper to convert the JP2s to GeoTIFFs and that works fine.
"geojasper -f in.jp2 -F out.tif -T tif"
However, when I use GDALWARP to change the EPSG:
"gdalwarp -overwrite in.tif out_a.tif -s_srs EPSG:3857 -t_srs EPSG:4326"
the resulting GeoTIFF's colors are washed-out when viewed in QGIS (still RGB, but very, very light, almost as if the brightness is set too high). Any ideas as to what might be the cause?
I will use a batch file to process all of the files once the problem is resolved, which is much, much faster than QGIS (surprised I can't select multiple files within QGIS to do this...).


texts: How to handle a dashed border line that becomes solid because it is shared

texts:
I am using Arc 10.3.1 and I am having a cartographic dilemma.  As shown in the picture below, I have dashed line representing borders of various city zones.  Where these border lines share borders with other polygons the resulting dashed line becomes much thicker and in some cases a completely solid line.  I can handle a slight variation in the size of the dashes but a solid line is unacceptable to me.

I have tinkered with changing the outline size and the outline dash type but each time when I export the lines are still merged. I have used definition queries to exclude the polygons in the middle, allowing the outer polygons lines to form the border for the interior ones, but there are several that I can not exclude because they create an outside border on the map that must be shown.

Is there a method to fix this little issue of mine?


texts: I want to rotate my North oriented point shapefile to 45%, is there any way to rotate it QGIS or MapInfo?


texts: How can I remove only small inner rings in PostGIS?

texts:
I am looking for a way to remove small inner rings in a large polygon.

I want to remove only the ones that are below a given surface and keep the bigger ones so that the the look (from a distance) is not affected too much.

If that could be done in a single SQL order, that would be great!


texts:
I am trying to set a categorised renderer for a point layer, according to its categories. I used this code:

vector = QgsVectorLayer("route to shape", "title", "ogr")
mi_rend = QgsCategorizedSymbolRendererV2()
mi_rend.setClassAttribute('TheAttribute') #Do I need one more attribute?
vector.setRendererV2(mi_rend)
QgsMapLayerRegistry.instance().addMapLayer(vector)

The problem is that the layer is loaded in the table of contents, nevertheless, it doesn't appear.
I have been looking for solutions in this forum but I haven't found any solution at the moment.


texts: I am trying to create a line style that consists of 2 lines. The first line (red) is actual border and the second line is a thick transparent line that should align it's bottom edge to the first line. When I zoom enough in QGIS, these lines render correctly, but when zoomed out, you can see that thick line is not rendered correctly. Problem is present even when I export map to PDF. Is there any solution for this?


texts: I am trying to create a code that will accept the input of a shapefile, and separate each of its features into separate shapefiles into a specified directory location. Can anyone help ? I am trying to use pyshp module, but it is only able to delete one file.


texts: How do I get the bounds of a layerGroup in Leaflet?

texts:
I have added multiple geoJSON layers to my map by first adding them to a LayerGroup and then adding the layerGroup to the map.

Now I want to call fitBounds on the map with the  combined bounds of all the layers, the problem is that layerGroup does not have the method getBounds

what can I do?


texts: WMS or ArcGIS MapServer draw geometry at point

texts:
I have a WMS service and use ArcGIS MapServer sometimes. These services connect to data source and get geometries as image.

I have a WMS or MapServer service that gets data from a database. But I want to show some temporary points on print view. So I need to show a point that my coordinates. I will send a point coordinate on WMS service or ArcGIS service it will show the point on response tile, but it will not be on database.

Is this possible?


texts: How to create a metadata json member for an openlayers layers object and display in extjs grid?

texts:
I want to do the following: add

metadata: "this is the layer metadata"

....to my openlayers layer config.

This should then be presented in an extjs grid. I can use

alert(layername.metadata);

to display the string inside the config member so I know it works.

To display this in the extjs grid I have used the following code

var myData = [];
var x=mapPanel.map.layers;
for (var i=0; i<x.length; i++) {
  // alert(mapPanel.map.layers[i].metadata);
    myData.push(mapPanel.map.layers[i].metadata);
}

var store = new Ext.data.ArrayStore({
        fields: [
           {name: 'metadata'}
        ]
        });

store.loadData(myData);

var myStore = new GeoExt.data.LayerStore({
    map: map,
    layers: mapPanel.layers
});

var grid = new Ext.grid.GridPanel({
        title: "Available Map Layers",
        region: "center",
        height: 300,
        viewConfig: {forceFit: true},
        store: store,
        columns: [
            {
                header: "metadata",
                dataIndex: "metadata",
                sortable: true
            }
        ]
    });

var metadataWin = new Ext.Window({
                title: "Map layers",
                layout: 'fit',
                width: '600',
                height: 'auto',
                border: false,
                closable: true,
                collapsible: true,
                x: 450,
                y: 100,
                resizable: true,
                closeAction: 'hide',
                plain: true,
                items: [grid]
           });
    metadataWin.show();

the grid show the correct amount of rows, but only the first letter is shown....this is probably a newbie javascript error.

Can anyone see where i´m going wrong?

yours,

Rob


texts:
I am trying to speed up a QGIS plugin I wrote by adding a spatial index to a vector layer. The aim of my plugin is to simplify a line network by merging some of its features, based on some criteria.

Short version: at each iteration, I delete two features in my layer and create a new one. This works fine. When I do the same thing in my spatial index, the changes do not seem to be taken into account.

My code:

Main function:

# Creating the spatial index (works well)
self.spatialIdx = QgsSpatialIndex()
map( self.spatialIdx.insertFeature, featureList )

Subroutine (feature1 and feature2 are to be deleted and replaced by mergedFeature):

# Update spatial index (doesn't work?)
self.spatialIdx.deleteFeature( feature1 )
self.spatialIdx.deleteFeature( feature2 )
self.spatialIdx.insertFeature( mergedFeature )

# Update layer (works well)
dataProvider.deleteFeatures( [feature1.id(), feature2.id()] )
dataProvider.addFeatures( [mergedFeature] )

When I try to use self.spatialIdx later to detect intersecting features (with self.spatialIdx.intersects( someFeature.geometry().boundingBox() )), I get the ids of allegedly deleted feature1 and feature2 (instead of the id of mergedFeature).

So, is there something obvious that I missed? Do I need to call a method to somehow "rebuild" the index? I also have a small doubt about the concept of global variable with self.…. Being new to Python, this approach may be wrong.


texts: How do I know if the Model Builder has updated.in ArcCatalog?

texts:
How do I know if the Modelbuilder has updated ?

So when I worked on Modelbuilder and then make some modified and added or subtracted tools. Sometimes you would forgot when you come back to it in a few days later and you have two same Model Builder and you don't know which one is the most recently ?


texts: rasterVis polygon under raster layer in R

texts:
I want to plot both a raster layer and a polygon layer using rasterVis::levelplot(), however, I want to have the polygon layer appear under the raster layer.

library(sp)
library(rgeos)
library(raster)
library(rasterVis)

r <- raster(xmn = 0, xmx = 1, ymn = 0, ymx = 1, nrows = 100, ncols = 100,
            crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
r[] <- rnorm(100 * 100)
r[sample(1:ncell(r), ncell(r) / 2)] <- NA

pt <- data.frame(x = 0.5, y = 0.5)
coordinates(pt) <- ~ x + y
proj4string(pt) <- projection(r)
pt_buff <- gBuffer(pt, width = 0.25)
levelplot(r, margin = FALSE) +
  layer(sp.polygons(pt_buff, fill = "grey40", col = NA))

I want the circle under the raster layer:


texts: Running Python script as input in QGIS 3.x Graphical Modeler

texts: I am building a model via QGIS 3.12 Graphical Modeler. There, I want to input Python code just like expressions, vector field,vector layer and use it as a parameter (or may be algorithm if it is possible to build). Any suggestions?


texts: How to determine what Level a new Extent would be, without zooming to it

texts:
In the ArcGIS Server JavaScript API 2.4, there are some related methods for getting and setting the map's extent.

map.getLevel() determine's the map's current level (when using a Tiled layer)
map.setExtent(extent) changes the map's extent to the new extent
map.setLevel(level) changes the map's level to the new level

Is it possible to determine what level the value of extent would equate to, before changing to that extent? That is, before I pass in map.setExtent(x1, y1, x2, y2) I want to know what level this would display, so that I can adjust the value as necessary.

I can set a listener on map.onExtentChange, verify the new level, and use map.setLevel(map.getLevel() + 1) to fix the extent if necessary - but this means the map may draw twice (once after map.setExtent then again after map.setLevel)

As far as I can see, there is no equivalent of extent.getLevel - is there a way to calculate this without actually changing the map level first?


texts:
We have published a geoprocessing tool that gives Raster dataset as output. we are using arcgis javascript api 4.0 version to visualize it on map. But when we use gp.getResultImageLayer , the event returns MapImageLayer. Now, we want to apply elevations to this MapImageLayer. This layer can be generated dynamically based on the input ascii file which we give input to the Geoprocessing tool.

Is it possible to apply elevation to MapImageLayer?


texts: Leaflet-routing-machine: Can't create route and/or destination marker from onClick on map

texts:
I have a map component created using react-leaflet, and have a marker function in it, which is supposed to add markers naturally on the map based on the latitude and longitude from the event object after the click.
  function addMarker(e){
    var coords = e.latlng;
    setMarkerPointsForRouting(markerPointsForRouting => [...markerPointsForRouting, coords]);
  }

Interestingly in the RoutingMachine component, I have a log which show the first click on the map i.e. the latitude and longitude. And I do get the first marker, but that's it. No end marker?
This is my map:
import React, { useState, useEffect } from 'react'
import { Map } from 'react-leaflet';
import LocateControl from '../LocateControl/LocateControl.jsx';
import MapboxLayer from '../MapboxLayer/MapboxLayer.jsx';
import Routing from '../RoutingMachine/RoutingMachine.jsx'

export default function MyMap({}) {
var [zoom, setZoom] = useState(8);
var [map, setMap] = useState(null);
var [animate, setAnimate] = useState(false);
var [markerPointsForRouting, setMarkerPointsForRouting] = useState([])
var [isMapInit, setIsMapInit] = useState(false)

  useEffect(() => {
  console.log("markerPointsForRouting ", markerPointsForRouting);
  },[markerPointsForRouting]);

  function addMarker(e){
    var coords = e.latlng;
    console.log("coords ", coords);
    setMarkerPointsForRouting(markerPointsForRouting => [...markerPointsForRouting, coords]);

  }

  function saveMap(map){
    if (markerPointsForRouting.length) {
      setIsMapInit(isMapInit => true)
    }
    setMap(map)
  }

  var MAPBOX_ACCESS_TOKEN = process.env.MAPBOX_ACCESS_TOKEN;

  var locateOptions = {
    position: 'topright',
    strings: {
      title: 'Show me where I am, yo!'
    },
    onActivate: () => {} // callback before engine starts retrieving locations
  };

  return (
  <Map animate={animate} zoom={zoom}  onClick={addMarker} ref={saveMap}>
    <MapboxLayer
      accessToken={MAPBOX_ACCESS_TOKEN}
      style="mapbox://styles/mapbox/streets-v9"
    />
    <LocateControl options={locateOptions} startDirectly />

    {isMapInit && <Routing map={map} latLng={markerPointsForRouting} />}
  </Map>
  )
}

And this is my Routing Machine component:
import React, { useState } from 'react'

import { MapLayer } from "react-leaflet";
import L from "leaflet";
import "leaflet-routing-machine";
import { withLeaflet } from "react-leaflet";

class Routing extends MapLayer {
  createLeafletElement() {
  const { map, latLng } = this.props;

  console.log("latLng ", latLng);
   var startIcon = new L.Icon({
    iconUrl: 'https://raw.githubusercontent.com/pointhi/leaflet-color-markers/master/img/marker-icon-2x-green.png',
    shadowUrl: 'https://cdnjs.cloudflare.com/ajax/libs/leaflet/0.7.7/images/marker-shadow.png',
    iconSize: [25, 41],
    iconAnchor: [12, 41],
    popupAnchor: [1, -34],
    shadowSize: [41, 41]
  });

  var endIcon = new L.Icon({
    iconUrl: 'https://raw.githubusercontent.com/pointhi/leaflet-color-markers/master/img/marker-icon-2x-red.png',
    shadowUrl: 'https://cdnjs.cloudflare.com/ajax/libs/leaflet/0.7.7/images/marker-shadow.png',
    iconSize: [25, 41],
    iconAnchor: [12, 41],
    popupAnchor: [1, -34],
    shadowSize: [41, 41]
  });

  var leafletElement = L.Routing.control({
    waypoints: null,
     createMarker: function(i, waypoints, n) {
       var marker_icon;
        if (i === 0) {
         marker_icon = startIcon;
        }
        else if (i == n - 1) {
         marker_icon = endIcon
        }
        var marker = L.marker(waypoints.latLng, {
         draggable: true,
         icon: marker_icon
        });
        return marker;
     }
  }).addTo(map.leafletElement);

  leafletElement.setWaypoints(latLng)

  return leafletElement.getPlan();
 }
}

export default withLeaflet(Routing);


texts: How to Add Vector Overlay to Clipped Image in Google Earth Engine?

texts:
I'm creating a clipped composite from MODIS Satellite's NDVI image collection. I need to add a vector overlay representing country boundary to each image in a collection.
My current failed attempts:

Creating a new image from custom image collection containing empty image with outline and satellite image. This resulted in outline only shown within clipped image (or only outline if outline image is first element in custom image collection array). The code is as follows:

var cntry_boundary = {...}; // feature collection
var geometry_to_clip_to = {...}; // a feature collection

var image_collection = ee.ImageCollection(...).map(function(image) {
          return image.clipToCollection(geometry_to_clip_to);
});

// Creating vector outline image
var empty = ee.Image().byte();
var outline_img = empty
  .paint({featureCollection: cntry_boundary, width: 1})
  .visualize({palette: '000000'});

var resulting_img = ee.ImageCollection([ image_collection.first(), outline_img ]).toBands(); // just to demo with one
Map.addLayer(resulting_img);

Other attempt was to try to apply mosaic function to custom image collection instead of toBands() to get combined image.

How can I create a new single image showing country vector outline as well as clipped image?


texts: Unable to activate cloned ArcGIS Pro default Python environment due to missing packages

texts:
To be brief, I am unable to activate cloned copies of the default ArcGIS Pro python environment, receiving an error message that required packages are missing.
I am using ArcGIS Pro 2.7.0

I am trying to create a new ArcGIS Pro python environment.
It seems the best way to do this is to clone the default environment, and then add packages using the GUI package manager or conda as necessary.
As seen in the picture below, I have 3 environments - sort of. The test environment was created by manually copy and pasting the default python environment folder to the root directory of my C:\ drive.
The arcgispro-py3-clone environment was created by clicking the clone default button at the top of the window seen in the image below.

Okay, so it seems like the environments were created correctly at first, but they're unable to be activated.
Lets try to activate the test environment. I click the radio button next to the test environment option.

Now I click OK in the bottom right.

A progress bar cycles for 5 to 10 seconds, and then I receive the following error: "Failed to set as active. Environment must have required packages for ArcGIS Pro."

I get the exact same error when trying to activate the arcgispro-py3-clone environment that was created using the Clone Default button as described above.
The error message makes me think that somehow my default arcgispro-py3 environment is missing packages. To try and fix that, I uninstalled ArcGIS Pro completely, as described in this ESRI blog post. Then I did a clean reinstall with all the default setup options, default install paths, etc.
I should also note I have tried restarting ArcGIS Pro multiple times, and have tried running as administrator as well.
ArcGIS Pro indicates that my default environment has 183 installed packages.

I have tried to find a list of the default packages required for ArcGIS Pro so that I could try and install everything to a cloned environment. I haven't been able to find a complete, or even partial list.
I have tried to google solutions for the exact error I'm getting, like ArcGIS Pro python environment "Failed to set as active. Environment must have required packages for ArcGIS Pro."
When I include quotation marks like that to search for my exact error message, i get 0 results on google.
I am at a loss at this point, but I think a first step in the right direction would be if someone could kindly examine their default ArcGIS Pro python environment, and verify if they have 183 installed packages like me. If someone was to report back that they had more than 183 packages, that would indicate that I do indeed somehow have missing packages.
Obviously, that wouldn't actually solve the problem. But it might kick-start the troubleshooting process for everyone...

I could post the output of conda list to provide a list of all the packages I have installed to the default python environment. However the output is many many lines, and I'm not sure what the preferred method of posting such a long block of text is here on StackExchange.